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PREFACE 


The Principal Investigator of this project is Professor Ivan I. 
Mueller, Department of Geodetic Science and Surveying, The Ohio State 
University. The Science Advisor is Dr. David E. Smith (Code 921, Geody- 
namics Branch), the Technical Officer is Mr. C. Stephanides (Code 942), 
and the Contracting Officer is Lillian E. Walker (Code 289), all of 
NASA/ Goddard Space Flight Center, Greenbelt, Maryland 29771. 

This report was originally presented as a dissertation to the 
Graduate School of The Ohio State University in partial fulfillment of 
the requirements for the PhD degree. 


ABSTRACT 


This investigation studies the possibility of improving the accu- 
racy of geodetic results by use of simultaneously observed ranges to 
Lageos, in a differencing mode, from pairs of stations. Despite their 
complexity and expensive computational requirements, orbital models are 
still not accurate enough to achieve geodetic parameter accuracies 
comparable to those of laser measurements. Parameters of interest here 
are the baseline lengths and the coordinates of the pole. Simulation 
tests show that model errors can be effectively minimized by simultaneous 
range-differencing (SRD) for a rather broad class of network-satellite 
pass configurations. To generate the required quasi-simul taneous range 
events, we compare the methods of least squares approximation with mono- 
mials and Chebyshev polynomials and the cubic spline interpolation. The 
latter seems preferable in the case of uniform and dense sets of data. 

Analysis of three types of orbital biases (radial, along- and 
across-track) shows that radial biases are the ones most efficiently min- 
imized in the SRC mode. The degree to which the other two can be 
minimized depends on the type of parameter‘s under estimation and the 
geometry of the problem. Sensitivity analyses of the SRD observation 
show that for baseline length estimations the most useful data are those 
collected in a direction parallel to the baseline and at a low elevation. 
Estimating individual baseline lengths with respect to an assumed but 
fixed orbit not only decreases the cost, but it further reduces the 
effects of model biases on the results as opposed to a network solution. 
Enforcing the simultaneity constraint ensures that the results are also 
free of possible biases introduced by inconsistent coordinate systems in 
which independent station determinations refer. Analogous results and 
conclusions are obtained for the estimates of the coordinates of the 
pole. 
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1. INTRODUCTION 


1.1 From Time Interval Measurement to Range Inference 

The laser systems which are deployed in satellite ranging are 
generally categorized as one type of electromagnetic distance measur- 
ing device (EDM). In a wider sense though, one can think of them as 
being one type of communications system. After all, in order to infer 
the distance defined by the instrument and the satellite, the trans- 
mitted signal must travel to the satellite and upon reflection return 
to the receiving system of the instrument. Assuming the velocity of 
the signal is known, the length of the signal path can then be computed, 
if the elapsed time is measured. The concept is remarkably simple; 
however, stringent accuracy requirements for the inferred distance 
demand that a great deal of calibration of internal errors and correct- 
ing for other systematic effects (e.g., atmosphere, retroref lector 
array biases, etc.) has to be done before the accuracy reaches accept- 
able levels. What this implies is that careful monitoring of the per- 
fonnance of every component of the instrument is required on a short 
as well as a long time basis, detailed recording of parameters descrip- 
tive of the environment, and, finally, tedious calculations to compute 
and apply the corrections. At times, much to the dismay of the scien- 
tist, even after all these corrections have been applied, the results 
do not reach the expected accuracy level. Human errors not excluded, 
the reason most of the time is the fact that the applied corrections 
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are in error themselves. Since none of the physical processes 
involved is known perfectly, we have to rely on laws and models based 
on observation in computing their effects on the observable. 

Naturally, the errors inherent in these models will propagate 
into our computations. In addition to this, the parameters needed to 
evaluate these models are also obtained from observation (e.g., tem- 
perature, pressure, clock offsets, clock rates, etc.) and therefore 
carry their own uncertainties due to observational errors. Even with- 
out having to go into a detailed description of the various steps 
involved in computing the final value for the station-satellite dis- 
tance, it should be by now obvious that an enormous number of factors 
have introduced uncertainties of different levels during this process. 
What is then provided to the analyst as an "observed range" is no more 
than a number cut of some computer software package along with a rather 
subjective estimate for its accuracy. 

1.2 Systems Accuracy; Some General Remarks 

It is common practice to implicitly assume that the error spaces 
of the various factors affecting a system are nearly orthogonal and 
that the interactions between them are therefore negligible. Although 
our present experience has not yet made a strong case against such 
practice, most statisticians [Scheffe, 1959] and metrologists [Eisen- 
hart, 1963] maintain that a system is best calibrated as a whole while 
in normal operation rather than in a piecemeal fashion or ideal condi- 
tions. We must clarify here that their definition of a system is the 
entity which consists of hardware and software components as well as 
operators and analysts. It is thus clear that quotations such as "our 
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present ranging capabilities to satellites" are rather meaningless 
unless they are supplemented with an explanation of the particular 
system we are referring to. Naturally, the satellite itself is a very 
important component, and even if everything else were to remain the 
same, changing the satellite can lead to order of magnitude changes in 
the system accuracy. 

1.3 Scope and Philosophy of the Investigation 

The scientist who is interested in a particular subset of param- 
eters will usually blame the inadequate hardware if the expected level 
of accuracy is not finally achieved. To a certain extent this may be 
true; however, hardware improvements are hard to come by and in most 
cases are very expensive. One then has to reassess the accuracy re- 
quirements and seek alternative solutions as well. The areas where we 
can look for improvements are the instrument design, the experiments or 
i.iission design, and the method of analysis of the collected data. 

This study addresses the third problem in connection with the 
estimation of interstation distances and variations in the coordinates 
of the pole from a geodetic satellite ranging system dedicated to geo- 
dynamics research: the Lageos system. As tar as instrument design 

improvements are concerned, it is unlikely that a geodesist can contrib- 
ute directly to any significant extent. It is important, however, that 
the geodesist has an understanding of a little more than its very basic 
principles since it is only thus that problem areas can be pointed out 
and subsequently looked at by the specialists. This knowledge will 
also help in communicating and exchanging information and ideas among 
the various science disciplines involved. If nothing else, it makes 
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it easier to foresee the miracles that the electronic gadgetry can 
or cannot achieve. 

The design of the experiment is an area where the geodesist will 
contribute more than anyone else involved. It is also the area where 
the most disappointment will be suffered. Even when the proposed 
design is a truly optimal one for the particular problem, what is final- 
ly implemented rarely bears more than a vague resemblance to what was 
originally conceived. Geographical, economical, political and other 
similar factors will usually reshape the design at the expense of op- 
timality. There is very little that one can do about this other than 
accept what is made available. 

Finally, the data analysis is the area where the geodesist can 
really experiment and innovate with practically no other than economi- 
cal limitations. The issue that this study addresses in that respect 
is whether a new method of analysis can be formulated which will mini- 
mize, if not completely eliminate, the effect the biases (inherent or 
introduced in the supplied data) have on the final results for the 
estimated parameters. The method proposed and investigated here is 
based on a liriear transformation of the range data to range-difference 
data. The range-differencing here refers to observations made simul- 
taneously from two ground stations to one satellite point as opposed to 
classical range-differencing of observations made from a single station 
to two consecutive satellite positions. To avoid confusion the pro- 
posed method will be hereinafter referred to as "simultaneous range- 
differencing" (SRD). 
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To summarize what has been said to this point, the scope of this 
investigation is the search for a new approach in analyzing range ob- 
servations which can be less expensive and time consuming while it will 
still provide estimates for the parameters that are of a higher or at 
least comparable level of accuracy to those obtained otherwise. The 
philosophy of the investigation can probably be best summarized in the 
cliche "The more expensive is not necessarily the better." 

1.4 Outline of the Proposed Simultaneous Laser Range-Difference 
(SRD) Mode 

Laser systems currently deployed in satellite tracking have 
recently been upgraded to accuracy levels where biases from systematic 
unmodeled effects prohibit us from extracting the full amount of infor- 
mation contained in the observations. Considering that the instrumen- 
tation quality improves at a faster pace compared to the available 
physical models, one can foresee that in the near future (when, for 
instance, NASA replaces all its lasers with third-generation models) 
the limiting factor for estimate accuracies will be the aforementioned 
biases. In light of these advances in the technological sector, it 
is only natural to look for new methods for the reduction of the obser- 
vations in ways that the effect of the biases can be kept well below 
the noi se 1 evel . 

The spectrum of geodetic satellite positioning techniques has 
been vastly enriched in the rather short quarter century life of this 
discipline. It has though remained polarized between two basic con- 
cepts: the geometrical positioning and the dynamical positioning. 

Over the years a number of hybrid techniques have also appeared to fill 
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the gap between the two dominant constituents of the spectrum. In the 
early days, simult'neous direction observations to satellites were suc- 
cessfully utilized to determine interstation directions [Veis, 1967; 
Aardoom et al., 1967; Schmid, 19741 using the formalism of geometric 
satellite geodesy and camera observations to balloon satellites such as 
Echo, Pageos, etc. Ranging systems never really participated in these 
solutions except for SECOR and C-band radar [Mueller et al., 1973] with 
rather poor quality observations. Apart from the fact that laser sys- 
tems were not as developed as they are today, the stringent requirement 
of the geometric solution for simultaneous observations from at least 
four and preferably more ground stations made the inclusion of laser 
ranging observations impossible. Still, geometric solutions are the 
only ones that do not rely on the dynamics of the satellite orbit and 
therefore the only ones which are not affected by their imperfections. 
Despite their large number of solve-for parameters and severe data dis- 
tribution requirement to avoid critical configurations and ill- 
conditioning [Blaha, 1971; Tsimis, 1973], their contribution to geodesy 
is vast and important. 

As models for the orbit dynamics improved, attention shifted 
rapidly to dynamic solutions and the geometric solution sustained a 
period of hibernation, only to be recently revived due to the develop- 
ment of the airborne laser ranging system by NASA [SGRS Workshop, 1979]. 
Full-fledged dynamic solutions are very expensive and involve thousands 
of unknowns with observations spanning several years and a number 
of satellites. Necessity, therefore, and the fact that the shorter 
the orbit the less time the orbital biases have to build up and corrupt 
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the solution, gave birl i to hybrid semi -dynamical solutions (short-arc, 
translocation, etc.). These techniques have been extensively applied 
in the case of Doppler satellite tracking observations which is natural 
if one takes into account how popular Doppler equipment has become in 
the scientific as well as the commercial sector [Brown, 1976; Kouba, 
1979]. 

In the case of laser observations, the equipment has been limited 
and mostly of observatory tyoe, difficult to relocate and for most uses 
too far apart to allow for significant amounts of simultaneously ob- 
served satellite passes. Furthermore, truly simultaneous events would 
be impossible to obtain since there is always an unknown synchronization 
difference between the clocks of the various stations. The launch of 
Lcgeos, however, has improved tremendously the geometry of the problem, 
since due to its high altitude (-5900 km) it is now possible that even 
intercontinental stations can co-observe th-’s satellite. As for the 
station mobility, the advent of TLRS I and the recent deployment of 
CLRS/TLRS II, as well as the international trend towards highly mobile 
and self-contained "observatories on whee I s/wings, " will soon allow 
for rapid deployment of instruments in almost any area that calls for 
it. Finally, when satellite time transfer becomes operational, it is 
hoped that global laser networks will be synchronized to better than 
100 ns compared to today's -1 us. Still, though, the use of laser 
ranging in a truly geometrical mode is highly unlikely due to its abso- 
lute dependence on weather conditions, a factor which is beyond our 
control (at least for the foreseeable future). That, however, has not 
stopped scientists from seeking alternate ways to improve the quality 
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of results obtained from the analysis of laser ranges to 
satel 1 ites. 

One of the most notable attempts to take advantage of geometry to 
improve the estimation of interstation distances from satellite laser 
ranging observations is summarized in [Latimer and Gaposchkin, 1977]. 
Their method, scalar translocation, takes advantage of coobserved satel- 
lite passes over the baseline under estimation. With rather poor data 
they have reported results that tend to be almost an order of magnitude 
better in accuracy and compare very favorably to independent estimates 
of the same baselines. Their success prompted us to undertake the 
investigation of using not only the coobserved part of the satellite 
pass, but, in addition, of converting the ranges to range-differences 
in hopes that they will be less affected by biases in the orbital model, 
the reference system, and in the observations themselves. Since there 
are not data taken specifically for this type of reduction technique, we 
had to select passes which had been coobserved by chance and then gener- 
ate simultaneous ranges from an interpolation of the recorded range ob- 
servations. Using then the generated simultaneous ranges from the end 
points of each station pair, we determine the simultaneous range differ- 
ence (see Fig. 1). These quasi-observables are then analyzed to obtain 
the minimum variance estimate of the baseline length. 

It is noteworthy that after this proposed investigation had been 
accepted by NASA, scientists at the Goddard Geodynamics Branch 
studied through an error analysis the advantages of using coobserved 
satellite passes in baseline determinations [Christodoulidis and Smith, 
1981], and their conclusions are in favor of this concept. Their 
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Data Collected at B 



Overlap Period for 


study indicates that even with the currently available instrumentation 
and physical models, a short arc of only five days with just six to 
eight passes coobserved by the baseline end stations can yield a 2o 
accuracy of better than one part in 2 x 10’ for baseline lengths ranging 
from 200 to 500 km. 

With such encouraging results when only implicit use of the geom- 
etry was made (only coobservation of the same portion of the orbit was 
required), it seems that our effort to make explicit use of it (simul- 
taneity--the same satellite position must be coobserved) is well justi- 
fied. Aside from this fundamental difference in the use of geometry, 
this investigation goes beyond baseline estimation and addresses the 
equally important subject of polar motion estimation from the same type 
of observables. 

As is well known and also discussed in this investigation, one can- 
not determine all station positions (and baselines therefrom) and the 
coordinates of the pole on the basis of the same satellite laser range 
data. To determine either of these types of parameters, the other should 
be known. The systematic errors in the latter do, of course, affect 
the estimation of the former. The choices here are rather limited. We 
could either fit some arbi trary functions to reduce the error growth-- 
udmittedly not a very attractive soluticn--or we can use the proposed 
simultaneous data reduction technique to minimize the effect of those 
unmodelled errors on the estimated parameters. 

The presentation of the material fcilows the natural sequence in 
which the questions arise and the results are presented, interpreted and 
discussed in the final chapters. The second chapter summarizes the 
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theoretical foundations required in this study. It discusses the refer- 
ence systems and frames as they are later used in the investigation and 
the orbital model for the satellite on which this study focuses most 
(Lageos). The third chapter is devoted to the observable, with a brief 
discussion on the instrumentation and an extensive investigation of tech 
niques for generating quasi-simultaneous range-differences when lacking 
simultaneous ranging events. The fourth chapter deals with the estima- 
tion process and related topics such as the estimability of the param- 
eters, the sensitivity of the observable in those parameters and the 
optimal network configurations for their estimation. 

The fifth chapter is a summary presentation of numerical experi- 
ments performed during the course of study of the proposed technique, 
using mainly simulated data and some real data which were available. 

The report of the investigation concludes with a chapter summarizing 
our conclusions and listing our recommendations for future research in 
this area. 
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2. THEORETICAL BACKGROUND 


2.1 Reference Systems and Frames 

One of the basic tasks in designing and carrying out an experi- 
ment is the collection of a set of rules which will properly describe 
the space in which the experiment takes place. We call this set of 
rules the reference system. 

( 

In almost all cases these rules are abstract in nature and pre- 
sent no means by which one could materialize the system in practice. 

That is, however, the ultimate goal, to be able to realize the concepts 
contained in the rules. The means by which we achieve this constitute 
the reference frame. 

The various refe*'P'^ce systems and the reference frames bear a 
relationship much like that between two volumes of an encyclopedia. 

Although the covers are the same, the content is quite different. In 
our case the same reference system can be realized in different ways, 
each one being a distinct reference frame. The distinction derives 
from thi different means by which the system is realized (stellar cata- 
log, planetary ephemerides, quasar source catalog, etc.) as well as 
from the different numerical values adopted in each case. Even though 
the concept of a reference system may imply some sort of uniqueness, 
the realization in practice may ha\e infinitely many variants. To 
remind ourselves of this arbitrariness we adopt the qualifier "conven- 
tional" reference systems and frames. We finally note that certain 
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assumptions and approximations made in the course of defining a sys- 
tem may render its validity as only an approximation of the ideal sys- 
tem and therefore in some cases the additional qualifier "quasi" needs 
to be used to indicate this. 

The nomenclature adopted here follows closely that proposed in 
[Kovalevsky and Mueller, 1981] since that is the most complete, self- 
consistent and internationally accepted one today. 

Although from the theoretical point of view one reference system 
would be sufficient for the description of the experiments, from a 
practical point of view we traditionally distinguish between space- 
fixed and earth-fixed systems. Both systems describe the same space 
continuum, though from a different point of view. The issue here is 
not the variety or number of systems to be involved, but rather the for- 
malism from which the rules defining these systems will be drawn. 

Should we follow the classical Newtonian formalism or should we follow 
Einstein's geometro-dynamics (general theory of relativity). 

Moritz [1979] suggests that a midcourse be sought. For all geo- 
detic purposes, the classical formalism amended with small corrections 
to account for relativistic effects can serve as a sufficient approxi- 
mation. It is acknowledged, however, that one place where the classi- 
cal formalism will fail is the fourth coordinate of the continuum: 
time. The reason is twofold: Einstein's curved space-time has defi- 

nitely discarded the (Newtonian) notion of a "universal" time, and, 
since time measurements today are the most accurate of all (less than 
one part in 10^^), there is no room for compromise or approximation. 
Each reference system has to have its own time scale associated with 
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it. This should come as no surprise, since as Minkowski [1908] him- 
self put it, "Henceforth space by itself, and time by itself, are 
doomed to fade away into mere shadows, and only a kind of union of the 
two will preserve an inc'^^^indent reality." For these reasons we have 
chosen to follow the theory of geometrodynamics in the operational 
definition of the required Conventional (quasi-) Inertial System (CIS). 
The corresponding Conventional Terrestrial System (CTS) is also defined 
in a consistent way, the former through the geodesic equations (or equa- 
tions of motion) of the satellite in the field generated by the surround- 
ing bodies (Sun, Earth, Moon), and the latter by an adopted set of dis- 
tances between a globally distributed set of stations--the CTS poly- 
hedron. Since working with these distances is quite cumbersome, a set 
of coordinates for these terrestrial stations is used, consistent with 
the distances within the error of measurement [Bock, 1983]. 

If our experiments were confined to events and processes well 
defined within the frame of the CTS, then we would not need to discuss 
it further. But we do refer to processes and events occurring outside 
the CTS, and we therefore need to ascertain that the above definition 
provides for accurate connection with other frames. With the origins 
of the two systems coinciding at the center of mass of the earth, the 
problem to be solved is that of the orientation of the axes of the 
crust-fixed CTS with respect to another frame which is observable with 
respect to the available observations. For reasons explained in 
[Mueller, 1980] and investigated in [Leick, 1978], we feel that the 
appropriate system is the Celestial Ephemeris System (CES). The third 
axis of this system is called the Celestial Ephemeris Pole or simply 
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the Celestial Pole (C). The defining principle for C is that it should 
have "no periodic diurnal motion relative to the crust (not the mantle!) 
or the CIS" [Mueller, 1980]. The adoption of this axis and the geocen- 
ter constrains the remaining two axes to a plane that is perpendicular 
to C and contains the geocenter. The motion of C with respect to the 
CTS third axis is described by the (observational ly determined) coordi- 
nates of the pole (Xp, yp) and the connection between the two is accom- 
plished through two well-established orthogonal rotations [Mueller, 
i969]. A remaining third rotation accounts for the earth's (i.e, poly- 
hedron's) sidereal rotation. This angle is modeled by a polynomial in 
time based on Simon Newcomb's expression for the right ascension of the 
fictitious mean Sun [Newcomb, 1898] and a linear-in-time component 
based on the mean sidereal spin rate of the earth. The true angle is 
obtained by correcting for the nutation in right ascension, known as 
the Equation of the Equinox (Eq. E) [Mueller, 1969]. A small correc- 
tion must also be added (determined observational ly again) to account 
for irregularities in the earth's spin rate. The resulting angle 
determines the instantaneous angular separation of two corresponding 
meridians of the CTS and the Celestial System, for example, the angle 
between their first axes. 

The relative motions of these frames are relatively well known 
in theory [Mueller, 1969]. The numerical models, though, are incom- 
plete and subject to continuous improvement as more observations be- 
come available and as their accuracies increase. In dealing with real 
data we have used the latest numerical models adopted by the lAU for 
use from 1984 onward [Kaplan, 1981]. In brief, the precession 
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formulation is that published in [Lieske, 1979], the nutation series 
are those derived by Wahr [1981], and the relationship defining the 
angular separation between the first axes of the CES and CIS (commonly 
known as the "Greenwich Apparent Sidereal Time") is taken from [Kaplan, 
1981]. The CIS-referenced orbits of the major perturbing bodies con- 
sidered in this investigation are taken from the JPL Development Ephem- 
eris 114 (DE114) and the corresponding Lunar Ephemeris 58 (LE58) 
[Standish, 1981, private communication]. 

The time scales used in this investigation are the Universal 
Time Coordinated (UTC) [Mueller, 1969] for data tagging purposes and the 
Barycentric Dynamical Time (TDB) [Kaplan, 1981] for the integration of 
the orbit. Both time scales are related to the International Atomic 
Time (TAI) scale, and therefore the relationship between the two is 
also known [Moyer, 1981a; 1981b j. 
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2.2 The Lageos Orbital Model 

The orbital model for a satellite consists of a set of formulas 
that determine the accelerations the satellite experiences within the 
frame of reference. As such, all the accelerations in our case will be 
referred to the CIS. For the particular case of Lageos, the accelera- 
tions which will be included are the terrestrial, lunar and solar point- 
mass effects, the terrestrial nonsphericity effects, the effects due to 
solid earth tides, the solar radiation pressure effects and an along- 
track acceleration of well-established magnitude but of as yet undeter- 
mined cause [Smith et al . , 1982]. Atmospheric drag effects are not in- 
cluded due to the high altitude of Lageos (-5900 km) and its very small 
cross-sectional area-to-mass ratio. All of the accelerations but the 
one due to the earth's nonsphericity are computed directly in the CIS 
frame. This deviation is due to the fact that the spherical harmonics 
that describe the earth's departures from a perfect point-mass body are 
given with respect to a body-fixed frame rather than an inertial one. 

To obtain an inertial expansion would mean to make those harmonics de- 
pend on time since they actually describe the anisotropic distribution 
of mass-density within the earth, which changes with time in an inertial 
frame due to the earth's rotational motions. For this particular accel- 
eration then, the determination results from a two-step procedure where 
first we determine the inertial acceleration in the CTS and then we ro- 
tate it into the CIS to make it consistent with the rest of the model. 
This particular acceleration is the only one involving the transforma- 
tion between the CIS and CTS frames. 
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The point-mass effects of the massive bodies in the solar system 
are determined partly on the basis of the geometrodynamical equations 
of motion, the Einstein-Infel d-Hoffman (EIH) formulation [Moyer, 1971], 
and partly on the basis of Newtonian theory. The former is used for the 
case of the solar, lunar and terrestrial effects, while the latter is 
reserved for the rest of the bodies in the solar system since their 
relativistic effects are too small to affect significantly Lageos' orbit 
[Moyer, 1982, private conmunication] . 

The determination of tidal accelerations is based on a series ex- 
pansion of the tidal ly induced potential by each of the considered 
celestial bodies (primarily the Moon and the Sun), truncated to the sec- 
ond degree. A more sophisticated model where the accelerations are com- 
puted on the basis of the tidal constituents (such as the ones derived 
by Doodson [1921]) is not required, since our goal here is not to study 
the tides per se, but rather to include their effects on the motion of 
the satellite in order to obtain a more realistic model for the actual 
motion. This formulation assumes a static earth and the same elastic 
response behavior of the earth for all orders within the same degree of 
the expansion. Although a dependence of this response on the frequency 
of the sustained tidal wave has been established [Lambeck, 1980; Gaposch- 
kin, 1981], for the purpose of the p>^esent study the simplified model 
is a sufficient approximation of the physical process. Since no special 
modeling of the ocean tides is included, it is justified to refer to our 
model as that for the solid-earth tides. Ocean tide effects are not 
modeled here since their effect is about an order of magnitude smaller 
than that of the solid-earth tides [Melchior, 1978]. As pointed out by 
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Musen [1973] though, both tides cause satellite perturbations with the 
same frequency spectrum and therefore aliasing is pssible. In fact, 
aliasing effects are responsible for the discrepancy between satellite- 
derived and ground-data-derived values of the effective Love number k 2 
that characterizes the earth's elastic response to the applied tidal 
attractions [Melchior, 1978]. For this exact reason the value k 2 = 0.27 
is deemed more appropriate for use here rather than that of k 2 = 0.30 
which is generally accepted as the more accurate estimate today [Lam- 
beck, 1980]. 

For the solar radiation pressure we have adopted a simple model 
based on the Sun's mean flux and a cylindrical shadow model as that given 
in [Cappellari et al., 1976]. A more sophisticated model would require 
tedious computations, and it is rather doubtful that the end result 
would justify the required effort. 

In the following sections we present the formulation for the com- 
putation of each individual acceleration and the corresponding variation- 
al equations. It should be mentioned here that the only unknov^ns in the 
orbital model are the position and velocity vectors of the satellite at 
the initial epoch of the integration. All other constants involved are 
assumed errorless although in practice one would solve for a number of 
them. For this reason, the variational equations which are presented 
here are only those that pertain to the unknowns; a more general set of 
equations can be found in [Cappellari et al., 1976]. The formulation is 
given in Cartesian rectangular coordinates since the special perturbation 
method of orbit determination has been chosen due to the complex accel- 
erations involved. 
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2.2.1 Orbit determination with the method of special perturba- 
tions. 

The determination of Lageos' orbit using the method of special 
perturbations [Cappellari et al., 1976] was dictated by the complexity 
of the orbital model and the stringent accuracy requirements. Although 
analytical methods can be more efficient and faster in orbit computa- 
tions, when the modeled physical phenomena become complex, those methods 
simply cannot handle them as accurately as numerical techniques and in 
some cases cannot incorporate them at all. In case of Lageos, no spe- 
cial formulation of the problem (time regularization, change of depen- 
dent variables, etc.) is required since the orbit is almost circular 
and very stable; we therefore chose the Cartesian formulation and coor- 
dinate time as the independent variable. 

A review of the relevant literature reveals that multistep numeri- 
cal integration methods are the most efficient and accurate for problems 
such as the determination of orbits. Furthermore, the integration of 
the original second-order differential equations is to be preferred to 
that of the reduced first-order system, since the higher-order set has 
a much larger region of stability. Higher accuracy can be achieved-- 
for a given stepsize--by increasing the order of the method; however, 
the higher the order, the less stable the process. We thus select an 
algorithm that permits variable order so that neither loss of accuracy 
nor instability can affect the solution. From the theoretical point of 
view the orbit of Lageos can be integrated with a fixed stepsize. Since, 
however, we use a variable order integrator, variable stepsize can result 
in some savings since as the order increases the stepsize is allowed to 
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increase too. A relative error test is used to determine whether to 
increase or decrease the stepsize and/or the order. The algorithm 
which fit our requirements and was conveniently available in computer 
coded form was that developed by Krogh [1969, 1970] of JPL. The struc- 
ture of this method consists of a modified Adams-Bashforth predictor 
and an Adams-Moul ton corrector of order one higher than the predictor. 
The mathematical formulation of these algorithms are given in [Krogh, 
1969], and the user required information in [Krogh, 1970]. 

Denoting by capital letters vectors defined in the CIS frame of 
the integration and by t coordinate time TDB, we can write symbolically 
the Lageos orbit equations of motion as: 


R 

where 


^PM ^ ^TD '^SR 



PM denotes point-mass effects 

N5 denotes nonsphericity of the earth effects 

TO denotes solid-earth tidal effects 

SR denotes solar radiation pressure effects 

AT denotes Lageos ad hoc along-track acceleration 


( 1 ) 


All of the above accelerations are functions of the satellite state 
vector and a number of model parameters (e.g., geopotential coefficients, 
etc.) which are assumed errorless in this study. Since the state vector 
at the initial epoch (initial conditions of the integration) is to be 
adjusted through the differential orbit correction (DOC) process to best 
fit the available observations, the variation in the initial conditions 
must be propagated to the epoch of observation state vector. This is 
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achieved by means of the transition matrix (or Jacobian of the state at 
the eooch with respect to the initial state) which is obtained through 
the integration of the variational equations of the state [Cappellari 
et al., 1976]: 


dt^VaP/ \3R/\3P/ \3^/dt\3P/ \3P/ 


( 2 ) 


where P denotes the vector of parameters of the DOC, in our case Rq and 
Ro, the initial epoch state vector, and the * on the last term indicates 
that the differentiation is carried explicitly with respect to P. Equa- 
tion (2) can be put in a more compact form if we define the following 
matrices: 



We can then write (2) as 

Y = A(t)Y + B(t)Y + C(t) 


(3) 


(4) 


Considering (1) we can write A(t) as 


A(t) = 


3R 


3R 


3R 


3R 


3R 


(5) 


where each term is a 3 x 3 matrix. Except for the relativistic part of 
Rpn^ and the R^j term in (1), no other acceleration depends on R. Since 
both these accelerations are extremely small, their contributions to 
the variational equations can be neglected with no loss of accuracy. 

The variational equations are quite insensitive to such small effects. 
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and in fact higner-order effects from the geopotential (for Lageos, 
above degree and order four) can also be neglected. This is important 
becai'-,*^ the computation of these effects is very tedious and computer 
time consuming. With no velocity-dependent terms in the model, we can 
set B(t) equal to the null matrix. In addition to this, C(t) can also 
be set equal to the null matrix because the accelerations do not depend 
explicitly on either RoorRo, which are the only elements of vector P. 
The variational equations that need to be integrated then take the sim- 
ple form 

Y = A(t) Y (6) 

with initial conditions (i.e., Y(t = 0)) 

Yo = [I I (7) 

I and ()) being the 3x3 identity and null matrix respectively. 

2.2.2 Gravitational acceleration from N point masses. 

In discussing the CIS system in Section 2.2, we pointed out that 
the chosen system will be realized in practice through the ephemerides 
of the bodies responsible for the gravitational field, as computed from 
the EIH formulation of their equations of motion. At this point, be- 
fore we give the explicit EIH equations, we must note again that this 
formulation is only accurate to an order 0 (1/c^), and it assumes no 
masses beyond our solar system. For the particular case of Lageos, we 
can further simplify the formulation by ignoring the relativistic ef- 
fects of all bodies except for the Earth, the Sun and the Moon. Fur- 
tnermore, since Lageos' mass is insignificant compared to that of the 
other bodies involved, we can safely assume that it does not contribute 
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to the generation of the field, and therefore it behaves as e true mass- 
less particle. Deviations of the assumptions on the internal structure 
of the perturbing bodies are not a problem either, since they can be 
taken into account through additional perturbing accelerations. The 
EIH equations have been derived by Moyer [1971] in a body-system- 
barycentric coordinate system for the general case of motion of body i 


in the environment of N bodies; 


uAr.-T.) 




2-1 kul 








E p .r . 
■ 2c^ j = l Ir.. 


where the required j-body accelerations are computed on a Newtonian 
basis from 


V. _ ^ (9; 

k^j 

Equations (8) and (9) are the basis for the derivation 0^ the 
equations of motion for Lageos. To obtain the acceleration of Lageos 


in the geocentered CIS frame, we use (8) to compute the Lageos acceler- 
ation with respect to the barycenter (>"pn/| in a second applica- 
tion, the acceleration of the geocenter with respect to the barycenter 
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“ B ^ I 

(Trm )q. The required acceleration for Lageos is obtained as 
their difference then: 



/— B 1 /— B» 

*'"PM *L ■ ^''PM 'l 


( 10 ) 


Since the first terms in (8) represent the Newtonian acceleration 
we can collect them together to produce a more illustrative form of the 
final equations. All of the other terms are divided by the speed of 
light squared, and we will use the notation (rpp|^ ) to indicate the 
sum of those terms, where the outer subscript indicates the event to 
which they refer, (L, Lageos; G, Geocenter). We can write (10) then as 





(r^-r) 



r 


“s 


(r^-r) 



+ 





-r) 



3 


'^IF-F 1 = 


i_Ml fJ_ ^ /- u, f- u. 
I- -13 ^ REL ^L ■ ^ REL ^G 

e‘ 




( 11 ) 


where the vector quantities on the right-hand side of the equation are 
all referred to the solar system barycentric coordinate system, and the 
differentiation is performed with respect to coordinate time, TDB. The 
above equation is written for the three major perturbing bodies, the 
Earth, the Sun, and the Moon as indicated by the subscripts E, S, and 
M respectively. Nonsubscripted quantities refer to the Lageos. This 
notation is illustrated in Fig. 2. 

For any other bodies of the solar system, only their Newtonian 
contributions are taken into consideration using the following expres- 
sion [Cappellari et al., 1976]: 
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2.2.3 Gravitational acceleration from terrestrial nonsphericity. 
The infinite spherical harmonics series which describes the geo- 
potential in space is given in [Heiskanen and Moritz, 1967]: 

00 n 

= — 1 [C 

r L m=0 

(15) 

It is implicitly assumed in writing (15) that the coefficients 
and are referenced to the (r,<}),A) coordinate system. The ac- 
celeration induced on the satellite by 4> is equal to its gradient, v^. 
Since the zero^^ term is considered in the point mass acceleration 


sin mA] P^^(sin<^)] 
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computations, we must subtract it from (15) which leads us to the "non- 
spherical" part of the perturbing potential V. The derivation of 
vV and subsequently is given in [Cappellari et al., 1976]. Because 
of the fact that the coefficients and are coordinate system 
dependent, the acceleration is obtained indirectly from the accel- 
eration of Lageos in the body-fixed system to which these coefficients 
refer. This computation involves the CIS to CTS transformation, and 
it assumes that these coefficients refer to the CTS. It is possible 
though that this may not be the case. If this is true then the classi- 
cal procedure must be revised and the harmonics must be transformed to 
the CTS. Formulas for such a transformation are given in [Kleusberg, 
1980]. Assuming that this transformation has been applied, and denot- 
ing by r^ the Lageos (inertial) acceleration in the CTS frame, then the 
acceleration in the CIS frame is 

= (SNP)^ r^ (16) 


J . 


where (SNP) is the CTS to CIS transformation [Mueller, 1980]. The 
acceleration vector r^ is obtained from vV as [Cappellari et al . , 1976] 


xz 

r^p 


_ 

r^p 

£_ 
L r2 


. ^ 
X 


_x 

r 

r 

1 

r J 


vV 


where 


(x,y,z) = (X,y,Z) (SNP) 


r^ = x^ + y^ T 


(17) 


(18) 

(19) 


and 
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( 20 ) 


The coordinates X,Y,Z are CiS referenced, and the x,y,z, CTS 
referenced. The 3x3 matrix in (17) t-'ansforms from the spherical coor- 
dinates r,(j),A to the Cartesian x,y,z coordinates. 

The associated va>^iational equations of state are given in 
[Cappellari et al., 1976] as 
■ - Far' 

^ = (SNP)^ (SNP) (21) 

.9RJ[^S 

and 

4 1 = [,] ( 22 ) 

3R 

with (p the 3x3 null matrix. 

Even if the coi.sistency of the coordinate systems has been 
assured though, there are still conditions under which the above equa- 
tions will give erroneous results. The reason is that certain low- 
degree harmonics, namely, C 20 . C 215 S 21 , C 22 » and S 22 are actually 

associated with the elements of the earth's inertia tensor [Heiskanen 
and Moritz, 1967; Nagel, 1976], and it can be shown [Nagel, 1976; 

Reigber, 1981] that tnrough these harmonics one can orient the refer- 
ence frame in which they are referred, with respect to the principal 
axes of inertia system as implied by these coefficients. The problem 
arises from the fact that since the earth is not a rigid body, its axis 
of figure (principal axis of maximum inertia) has body-fixed motions 
due to seasonal mass redistributions (free motion) and a diurnal motion 
due to the tidal bulge (forced motion). The free motion is rather small, 

t 

9 


29 


with an amplitude of at most 2m and a period close to the Chandlerian 
('430 days), but the forced diurnal motion can reach amplitudes of 
±60 m! These facts have been known for some time [McClure, 1973; 

Nagel, 1976; Leick, 1978; Moritz, 1979], but a combination of improved 
measuring accuracies, more stringent requirements for accurate results, 
and more dense (frequent) observational records have almost reached 
the point that one cannot afford to continue ignoring them for much 
longer [Tapley, 1982], 

In the present study we have adopted a time variant nature for 
the C 21 and S 21 coefficients, computing their values at each epoch on 
the basis of the adopted C 20 value and the coordinates of the pole, 
using the formulation of Reigber [1981], 

2.2.4 Solar radiation pressure acceleration. 

Photons emitted from a radiating source at a frequency v possess 
an energy hv, where h is Planck's constant. Since they travel in space 
with the speed of light, their momentum can be computed as the ratio 
of their energy to their speed. When a massive body travels through 
their continuum, the impinging photons transfer part of their momentum 
to the body upon impact. It is obviously impossible to know the exact 
number of photons that collide with a satellite, but we can determine 
the total force if we know the distribution of photons in the continuum, 
i.e., their flux. This flux, however, varies, and its instantaneous 
value can only be obtained through observation. For the sun, for 
instance, solar flares can significantly increase the flux, hence the 
photon count. Nevertheless, for the present study adopting a mean 
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flux is sufficient, since the effect on Lageos is quite small already. 

To determine the acceleration of the satellite due to solar radi- 
ation pressure, we first determine the force exerted on the satellite. 
This force is proportional to the satellite's area of impact A, the 
solar flux S, and a radiation pressure coefficient C^, and inversely 
proportional to the speed of light c, and the sun-satellite distance 
|R^-R| squared. The radiation pressure factor Cp is dependent on the 
reflectance characteristics of the impact area. This radiation pressure 
coefficient has been sufficiently accurately determined for Lageos at 
1.1729. The magnitude of the resulting acceleration on Lageos can be 
written, then, as [Cappellari et al., 1976] 


£ = 1 ^R ^ 

m c |Rc-Rl^ m 


(23) 


where S is the mean solar flux at one astronomical unit (1 AU), and F 
is the force exerted on the spacecraft. For convenience S/c is denoted 
by P^, the solar radiation pressure constant, and it gives the force 
exerted on a perfectly absorbing body (C[^=l) at a distance of 1 AU 
from the sun. We have adopted the value of about 4.626 x 10"^ Newtons/ 
m^ for 

Since this acceleration is the result of a collision, its direc- 
tion is opposite to that of the impinging photons, i.e., in the direc- 
tion of R-R^. We can then express the acceleration in the CIS coordi- 
nates as 



( 24 ) 
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Because the satellite is not always exposed to solar radiation, its 
motion can eventually bring it within the earth's shadow; equation (24) 
must be applied only part of the time in the computation of an orbit. 

To check whether the satellite is in or out of this shadow, one can go 
to various degrees of sophistication, from the simple cylindrical 
shadow model to the complex conical one, discriminating between umbra 
and penumbra regions [Baker, 1967]. We adopt the sample cylindrical 
model and we introduce the eclipse factor y with only two possible 
values: 

Y = 1 the satellite is in sunlight 

> = 0 the satellite is in the shadow 

To determine the value of y» the following computational checks 
are performed at each step: 

(a) If R-Rg > 0, then y = 1 |r x R | 

(b) If R*R^ < 0, then compute j-^ | ~ 

IR X Rj-l 

If D = ar >0, then y = 1 

l^sl ■ 

|R X Rjl 

If D = = — a_ < 0, then v = 0 

iRsI ' 

The logic behind this check is best explained in Fig. 3. 

If one wants to be absolutely rigorous in the computation of R^p, 

then the earth is not the only body to consider. The lunar shadow 

should also be considered. Furthermore, the above formulation considers 
only the effects of direct solar radiation. If radiation pressure is a 
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Fig. 3 Cylindrical shadow geometry. 


major perturbation in the motion of a satellite (unlike Lageos, A/m 
'0.007), then in addition to direct effects, one must consider the ef- 
fects of diffusely and specularly reflected visible radiation from the 
earth and also the opposing effect of promptly and delayed emitted 
infrared radiation from the body of the satellite itself [Baker, 1967]. 

The associated partial derivative of (24) that contributes to the 
variational equations of the state are given in [Cappellari et al . , 
1976] as 


. 9R . 
and 

. 9 ^ . 



(R-R3)(R-R3)^ 





where I and 41 denote the identity and null matrix respectively. 


(25) 
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2.2.5 Tidal acceleration. 

The attraction of the Sun and the Moon on the deformable Earth 
creates a time-varying (due to the relative motion of the bodies) poten- 
tial field which in turn induces additional accelerations on an Earth- 
orbiting satellite. This tidal potential can be expanded on the Earth's 
surface as an infinite series, as we have already discussed, though only 
the second-degree term is of any significance in this study. Accord- 
ing to A.E.H. Love [1911], the earth not being a perfectly elastic body, 
the actual response of the earth can be obtained from this potential 
as a fraction of it. The ratio of the induced potential to the tidal 
potential was designated by Love as k, a number whose value depends on 
the elastic properties of the earth. Further dependence of k on the 
frequency of the tidal waves has been established today, but for this 
study we will simply use the effective Love number k 2 with a numerical 
value of 0.2748, consistent with satellite results when no modeling of 
the ocean tides is included. 

The tidal potential on the earth's surface can be written then 
in the CIS frame as follows 

Ut = — ^ P 2 (cose) (27) 

where R. is the distance between the centers of mass of the earth and 

D 

the disturbing body b, is a mean radius for the earth, P 2 is the 
Legendre polynomial of degree 2, and e is the angle defined by the direc- 
tion of the geocentric radius of the evaluation point and the geocentric 
direction to the disturbing body. On the basis of (27) and using 
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Dirichlet's principle, we can obtaifi the potential that perturbs 
the satellite's motion as 


2 



where we have used the inner product between vectors to give a more 
manageable expression for Palcose). From (28) and the fact that the 
tidal acceleration on the satellite is equal to '?Uq^(R), we obtain 




5(u^-u)2]u + 2(u^-u)ujj 


(29) 


where u^ and u are unit vectors in the direction of the tide producing 
body b and the satellite respectively. 

From (29) we can obtain the contribution of this acceleration to 
the variational equations as 


.9R Jl 


f k2 


Pb 


|R 


[35(u^-u)2- 5]u u^ + 2 u^ uJ 


+ [1 - 5(u^-u)^]I - 10(u^^-u)[u uj + Ujj u^] 


(30) 


and 


3R 


TD 


dR 


[*] 


(31) 


where I and (}> are the 3x3 identity and null matrix respectively. The 
derivation of this equation can be found in Appendix A. 

The complete effects on the satellite are obtained through a 
summation of each body's effects over the number of bodies considered. 
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2.2.6 Lageos empirical along-track acceleration. 

As more dense and accurate observations became available from 
Lageos, its orbit was determined with increasing accuracy to the extent 
that an unmodeled secular acceleration causing a mere 1.1 mm/day de- 
crease of the semi-major axis could definitely be identified in all 
orbital fits [Smith and Dunn, 1980]. The mysterious acceleration, al- 
though small, alarmed scientists since its definitive existence jeopar- 
dized the long-term stability of Lageos' orbit and scientific interest 
was stirred up to attempt to give an explanation as to what is the phys- 
ical cause for it. In most cases, the research has been concluded with 
a rejection of the assumed physical phenomenon as the possible cause, 
which in science is sometimes a more valuable contribution thon finding 
the cause itself. Two most recent publications give an example of the 
diverse directions in which a solution has been sought. 

Rubincam [1980] has examined a number of physical processes char- 
acteristic of the upper atmosphere and others such as gravitational 
resonance, gravitational radiation and terrestrial radiation, only to 
conclude that of all processes examined only drag due to charged and 
most likely neutral particles can be a plausible cause for the acceler- 
ation. If the relevant theories were further developed, then a more 
definitive answer could be given. For the time being though, his main 
argument for accepting this explanantion is that it also helps solve 
the so-called "helium problem" [Chamberlain, 1978]. 

Szebehely [1981] on the other hand, keeping in line with his in- 
terest in the inverse problem of celestial mechanics, attempts to find 
the potential that would generate the unexplained acceleration. On the 
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basis of a circular orbit and two-body motion he shows that this poten- 
tial is proportional to the spiral function of Lageos (i.e., the equa- 
tion described by Lageos on its orbital plane) and inversely propor- 
tional to the square of its geocentric radius. One could then con- 
ceivably examine processes that are candidates for the solution of the 
problem and identify which of them could give rise to such a potential. 

Since there is no definitive answer to the problem yet, we chose 
to model the unexplained acceleration in the same empirical way that 
the GSFC researchers do. From the full amount of available observations 
on Lageos, a set of monthly values of the magnitude of this acceleration 
has been compiled. We have adopted as definitive values those shown in 
Table 1, communicated to us by Mr. Ron Kolenkiewicz of GSFC, covering 
the period of time that was of interest to us. In a recent report 
[Dunn, 1982], the average value of this magnitude over the first five 
years of Lageos data has been reported to be 2.86 x 10 m/s^. This 

value can be used as a standard of comparison for the monthly variations 
of the effect. 


The empirical model assumes that a restraining force is acting 
against the motion of the satellite, giving rise to an along-track (in 
the direction opposite to that of the velocity vector) acceleration 



One then can express this acceleration as 




(32) 


where a is the magnitude of the acceleration taken from Table 1 for the 


relevant time period. No contribution to the variation of the state is 




Table 1. Lageos Along-Track Acceleration Magnitude 


Time Period 

Acceleration Coe 

October 1979 

4.3 

Novembe*' 

4.4 

December 

3.2 

January 1980 

1.6 

February 

0.4 

March 

0.8 

April 

3.8 

May 

4.3 

June 

3.2 

July 

3.5 

August 

3.4 

September 

3.3 


computed due to the fact the effect is too small to bring any signifi- 
cant change in the variational partial derivatives. 
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3. THE OBSERVABLE 


3.1 Hardware Components 

The operational principle of a laser ranging system was already 
mentioned in the introduction of this study. The range is inferred 
from the round-trip flight time of a pulse emitted from the station, 
reflected on the satellite and received back at the station. It is 
thus obvious that the system consists of two major components, the 
ground-based active instrumentation and the spaceborne satellite target. 

3.1.1 Laser ranging instrumentation. 

The basic subcomponents of the ground-based part of the SLR sys- 
tem are the laser cavity, the timing equipment, the detection equipment, 
and the pointing system. Other secondary components are mini -computers , 
ancillary data collecting devices, calibration instruments, etc. 

We will restrict ourselves here tc jescribing the present capabil- 
ities of the available systems in terms of the observed ranges' accura- 
cy and pointing out the sources of systematic errors affecting them. 

For a more detailed explanation of these errors consult [NASA, 1980]. 

From the laser cavity itself, the major source of error is the 
distortion of the wavefront. It can map into the range as high as 6.0 
cm, but for the last (third) generation systems is no larger than 0.5 cm. 
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The receiver system contributes a much higher error which is a result 
of the photodetector, the discriminators and the delays from the various 
cables involved. Depending on the quality of the system, it varies 
between 1.5 - 8.0 cm. The timing system introduces systematic errors 
in the data for two reasons. The time interval counter frequency drift 
affects the measurement directly. The frequency standard that keeps 
absolute time and is used in determining the epoch of the measurement is 
also affected by similar factors, and it indirectly degrades the accu- 
racy of the data by supplying biased epochs. The total effect is no 
larger than 1.0 cm though. 

Finally, the calibration process itself, although by definition 
its purpose is to remove such errors from the measurements, can be in 
error too. In the older systems that calibration was performed through 
test measurements before and after a ranging session; the contribution 
to the error budget could be as high as 1.0 cm. New models which are 
capable of doing a real-time calibration have almost eliminated this 
error source [Silverberg, 1981a]. 

Although tiie above errors are systematic for a particular station, 
we cannot claim that they are systematic among stations. In fact it is 
more likely that they would be random in that sense. It may be that in 
some cases simultaneous range-differencing between two stations will 
result in a minimization of the above errors, but we cannot safely 
assume that this will be the case always. 



3.1.2 Laser Geodynamic Satellite (Lageos). 

The spaceborne segment of a satellite laser ranging system con- 
sists of a satellite which is equipped with corner cube reflectors (CCR) 
which reflect light at the incident direction. One such satellite which 
is dedicated to laser ranging is the Laser Geodynamic Satellite-- 
Lageos (7603901). We describe this satellite only because as its name 
hints it has been launched to support geodynamics research and therefore 
it is the optimal target for such purposes. 

The utility and merits of a target satellite such as Lageos are 
best summarized in [Johnson et al., 1976], the abstract of which we 
quote: 


The fundamental concept of Lageos is a long-lived, dense, 
electrically and mechanically inert spherical satellite with its 
surface speckled with retroreflecting cube corners, designed such 
that range measurements between duly equipped laser ground sta- 
tions and the satellite are possible with an ultimate accuracy of 
2 cm when data from a single satellite pass are appropriately 
averaged. The Lageos concept requires that the satellite be 
placed in an orbit for which an ephemeris can be determined 
ultimately to 5 cm rms uncertainty for a 24-hour arc. These 
required satellite characteristics should allow the several 
geodynamic motions experienced by ground stations to be deter- 
mined typically with 2 cm accuracy. 

Lageos is a sphere of 0.59988 m diameter, with a mass of 407.821 
kg, which results in a mass-to-area ratio of 1.44 x 10^ kg/m^ and 
therefore small sensitivity in solar radiation pressure and aerodynamic 
drag forces. The outer shell of the sphere is made of aluminum with a 
matte exterior surface which results in Lageos appearing as a star of 
visual magnitude 12-13. The sphere contains a core drum made of beryl- 
lium copper to provide the desired weight. 
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The surface of the satellite is speckled with corner cube retro- 
reflectors arranged as uniformly as possible in eighteen "rows" or 
"rings" and two single CCP's at each end of the satellite's axis of 
symmetry. There is a total of 426 high quality, precision manufactured 
and tested OCR's [Fitzmaurice et al . , 1977]. 

Following the basic concept of the Lageos system, its orbit was 
selected in such a way that perturbations caused by solar radiation 
pressure, drag, and high frequency poorly determined gravitational con- 
tributions were minimized. Nominal characteristics of the orbit, as 
determined on the basis of initial acquisition data obtained with the 
Baker-Nunn camera and laser system of SAG are given in Table 2 [Pearl- 
man et al., 1976]. Lageos has ever since been consistently tracked 
with the international laser network maintained by NASA, SAG, and indi- 
vidual country agencies. Its ephemeris quality has improved consider- 
ably as improved ground instrumentation has been deployed and as the 
parameters for the force model have been obtained with similarly better 
accuracy. It is claimed [Lerch and Klosko, 1982] that we can now model 


Table 2 Lageos Nominal Grbital Elements 


Epoch 

June 7.G, 1976 

Incl ination 

109?8585 

Eccentricity 

G.GG3929 

Apogee Height 

5941.9 km 

Perigee Height 

5845.4 km 

Period 

225.47G6 min 

Semi major Axis 

12 271. 79G km 
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the motion of Lageos over a month's time to better than 0.5 m. Never- 
theless, it is natural to seek orbital models that have accuracies com- 
parable to that of our observations which to date have reached the im- 
pressive 2-3 cm level. Once we achieve parity in that respect and sub- 
sequently maintain it, we can say that we have exploited the full poten- 
tial of the Lageos system. 

3.2 Generation of Simultaneous Range-Differences 

The most tedious part of this investigation, aside from the writ- 
ing, debugging and testing of the required orbital program, was the part 
pertaining to the generation of SRD's from the voluminous amount of 
laser ranging data provided by NASA/GSFC. The above statement though 
requires clarification, since the normally required effort is much less 
than what one might surmise from it. Most of the work involved can be 
characterized as "hunting" for suitable data. The fact that no specific 
campaign was ever devoted to the application of this new technique 
made data collected from stations coobserving the same pass hard to come 
by and only by coincidence. It must be noted here that this task was 
enormously facilitated by computer software kindly made available to us 
by Mr. Ron Kolenkiewicz of NASA/GSFC [Kolenkiewicz , 1980, private com- 
munication]. 

In the following sections we summarize the procedure that has been 
followed herein for the generation of the SRD's. We initially describe 
the range data set as originally obtained, the corrections that w'ere 
already applied to the ranges and corrections which we had to apply in 
addition. For the sake of completeness, we describe the details of the 
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selection process, although this in effect amounts to an explanation of 
the NASA-supplied software. Finally, since we would eventually need to 
interpolate the selected observations in order to generate the SRD's, 

W3 have investigated three popular approximation techniques, truncated 
power series, Chebyshev polynomials, and cubic splines, in order to 
determine which would suit our purpose best. This by no means should 
be taken as an exhaustive treatment of the problem since it is a re- 
search topic in its own right and one can probably spend several years 
testing various approximation methods each with its own merits and draw- 
back? as wel 1 . 

An investigation into the required relativistic corrections in 
the measured ranges and the corresponding SRD's showed that this correc- 
tion is for all practical purposes negligible for the latter. The for- 
mulation and derivation of these corrections along with the station 
dependent tidal corrections are given in Appendix B. 

3.2.1 Description of the original set of range data. 

Prior to the discussion of this section's topic, some clarification 
on terminology is in order, as far as the types of data encountered in 
laser ranging is concerned. We will refer to the data collected at the 
individual laser stations as "raw data." Besides the observation itself, 
each raw data set contains a great amount of ancillary data that is re- 
quired in the preprocessing of each observation. Raw data are not use- 
ful to an analyst unless the person is familiar with the operational 
details of the NASA- or SAO-supported networks of stations. Instead, we 
obtain the "preprocessed data" which comprise the actual observations 
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with certain corrections applied, and with other corrections having 
been computed but not applied. The latter are included in the dissemi- 
nated data set with appropriate indicators that inform the analyst of 
what has been included already in terms of corrections and what yet re- 
mains to be applied. Such practice is particularly helpful in case 
someone would prefer to compute the corrections anew with different 
models or revised constants. The quality of the raw data depends heavi- 
ly on the quality of the laser and timing instrumentation at each sta- 
tion. The preprocessed data, though, go through a number of processing 
stages and their quality is determined not only by that of the input 
raw data, but additionally by the integrity of those processes which 
they undergo. It is rather unfortunate that no scrupulous account has 
yet been published for either of these two stages of data handling. 

For both, the initial in situ processing of the collected data, and for 
the subsequent further refinement at the central computational center 
of NASA at &SFC, the documentation is outlined in a single short docu- 
ment [Carpenter, 1978] distributed to project investigators or obtained 
on request. 

The gist of this document is that there are in general four cor- 
rections which need be applied (and are) to the observed ranges; fixed 
threshold to peak (return) signal offset, instrumental calibracion cor- 
rection, atmospheric refraction correction, and satellite center-of- 
mass offset correction. Some other corrections might be done addition- 
ally to compensate for peculiarities of individual instruments (e.g., 
mount axis offset correction for X-Y type mounts). The first correction 
is obtained from the output of the waveform digitizer and the registered 
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travol time of the pulse. The internal calibration correction is deter- 
mined from ground-target ranging results collected just prior and imme- 
diately after a satellite ranging session. For the latest generation 
of laser instrumentation these two corrections are determined and ap- 
plied in real time [Silverberg and Malevich, 1978; Silverberg, 1981a]. 

The atmospheric refraction correction is computed from a formula given 
in [Marini and Murray, 1973] using the pertinent information collected 
at the station during the tracking session and transmitted along with 
the raw data to the computing center. The center of mass offset correc- 
tion is determined for each satellite before its launch. It is meant 
to refer the observed range to the point which is considered the center 
of mass of the vehicle and to which the computed orbit state vectors 
refer. For spherically symmetric satellites, such as Lageos, the correc- 
tion is a constant (0.24 m for Lageos), but for satellites with compli- 
cated figures (e.g., those which carry solar paddles) and whose center 
of mass location depends on the relative position of moving parts (such 
as the paddles), the computation of the correction depends on several 
factors and sometimes it can only be approximately known. This is one 
more reason why satellites such as Lageos are ideal for the precise 
positioning required in geodynamics research. 

When all the corrections have been applied, the data collected 
over each satellite pass are fitted to a polynomial to identify blunders 
and edit suspect observations. Nominally, a fifteenth-degree polynomial 
is used; however, for passes where less than 60 ranges are accepted, the 
degree of the polynomial is reduced to one-fourth the number of accept- 
able data points. The editing is performed by comparing the difference 
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between the observed range and that computed from the fit to the 
rms deviation of the fit. Data points which differ by rnore than three 
times this standard deviation are edited and the process repeated. 
Carpenter [1978] assures us that for relatively wel 1 -di stributed data 
points, this process converges quickly and indeed eaits only spurious 
data. It is not the goal of this investigation to examine the correct- 
ness of these subtasks, but it is worth noting that the so-called 3o 
rule has been contested on various occasions [Uotila,1973 and 1976], and 
the implicit assumption that the analyzed data errors are normally dis- 
tributed (for the rule to have any meaning at all) is made rather out of 
convenience than on the basis of some theoretical justification. At any 
rate, discussions at the meetings of the parties involved indicate an 
awareness of the problems in data preprocessing, ana it is safe to say 
that major revisions can be expected in the near future in several as- 
pects of this task. 

The corrected ranges and the applied corrections are archived at 
the National Space Science Data Center (NSSDC) from which copies may be 
obtained upon request. It is this type of data that we have obtained 
for the numerical tests of this investigation. The format in which the 
data have been encodrid on the magnetic tapes for transmittal is commonly 
known as the "Lageos binary format" [Putney, 1980], and it replaces the 
previously used (similar) Geos format. 
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3.2.2 Selection of ranging data on simultaneously observed 
satellite passes. 

The distributed preprocessed ranging data are obtained on magnetic 
tapes, arranged in files, one for each month of the period covered. 
Within each file, the data are arranged in chronological order. For 
each set of monthly data, a catalog is constructed showing the stations 
which have collected observations over that month and the number of ob- 
servations collected. The catalog cor.. dins also a pass-by-pass break- 
down of the data and the beginning and ending epochs for each station 
having observed a certain pass. On the basis of this catalog and the 
known geographical location of the ooserving stations, the station pairs 
which are likely to have sufficient numbers of observations on the same 
portion of a satellite pass are determined. 

If the number of coobserved satellite passes and the distribution 
of the observations seem promising, the next step is the actual deter- 
mination of the overlapping observational periods and the number of ob- 
servations collected by each station of the pair over those periods. 

This is accomplished by examining the data with the OVERLAP software. 

Due to either equipment or data failures, occasionally a pass is inter- 
rupted by gaps due to the missing or edited data. These gaps result in 
an uneven data distribution over the pass and can cause problems at a 
later stage, when simultaneous ranges are to be interpolated. To allev- 
iate this problem, the software checks the duration of these gaps in the 
selected passes, and if they are larger than an allowed value, the pass 
is broken down into subsets of data. A 30-second maximum gap has proven 
reasonable as it does not cause problems in the interpolation, and at 
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the same time it does not result in extremely small subsets that would 
be hard to interpolate due to insufficient data points. At this point 
the results of the overlap analysis are examined as to the data content 
and distribution; periods with less than ten data points are rejected 
and so are the data from station pairs that have no significant total 
amount of observations over the whole month. 

The remaining identified periods are now used to select the actual 
observations out of the original data set. At this stage, the distribu- 
tion of the selected data within the overlap period for each of the 
paired stations is determined. In almost all cases, and for obvious 
reasons, the station with the largest nuinber of observations also has 
the best distribution. It is important to known this detail, because 
only for one of the stations need the ranges be interpolated at the 
epochs that the alternate station has observed. We therefore choose to 
interpolate the ranges from the station with the best data distribution 
in order to keep approximation errors as small as possible. The end 
product of the selection process is two files, one containing the select- 
ed observations, and one containing a data-station directory. Included 
in the directory are the endpoint epochs fo>^ each batch of data consti- 
tuting a pass or a portion of it, the identification numbers of the sta- 
tions coobserving the pass, and an indicator that determines for which 
of the stations the ranges will be interpolated and for which the actual 
observations will be used. 


49 



3.2.3 Functional representation of the data selected for 
SRD generation. 

We discuss here three different approaches for the approximation 
of the station-satellite distance at a number of given epochs, on the 
basis of observations of this distance at other instances spanning the 
approximated interval. One restriction in our choice from available 
methods is the arbitrary distribution of the observed ranges over the 
satellite pass. Having no control at all over this factor, we have but 
to eliminate outright the possibility of using some very simple, effi- 
cient and accurate methods such as the Chebyshev approximation [Dahl- 
quist and Bjbrk, 1974] or, for that matter, any other method which re- 
quires that the given base points correspond to a particular set of 
values of the free variable. Of tne remaining viable methods, we have 
chosen to study two which either by virtue of their simplicity or their 
accuracy properties have attained widespread popularity among those who 
analyze experimental measurements. These methods are the least squares 
approximation and the cubic spline interpolation. There is a vast lit- 
erature for both topics; one, however, is almost obliged to refer to the 
classical text of Davis [1975] for the first, and the concise but practi- 
cal treatise of Spath [1974] for the second. As far as the least squares 
method is concerned, we have investigated the appl ication of the method 
with two different sets of base functions: monomials and Chebyshev 

polynomials. To test the three methods in the absolute as well as the 
relative sense, the following experiment was performed. 

We chose two sets of ranges, each set distributed over an interval 
of about 15 minutes, typical of the overlapping periods encountered in 
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this study. One set is chosen to have a high data density, with observa- 
tions being made at a rate of about 1 pps. The other set is one of very 
sparse data, with observations about 20 seconds apart on the average. 

The sparse data set was collected at the Orroral (7943) SLR station in 
Australia, while the dense one comes from the Goldstone (7115) SLR sta- 
tion in California. Figs. 4 and 5 are the plots of these ranges versus 
time. The vertical bars indicate the actual data points. From these 
two sets of ranges, we selected a number of data points evenly distrib- 
uted over the entire pass to form a "ground truth" data set. These 
selected ranges were eliminated from the original data sets. The re- 
maining observations were subsequently used in the approximation process, 
using all three different methods. The results are then used to approx- 
imate the station-satellite distance at the "ground truth" points pre- 
viously selected. Comparing the true values with the interpolated ones 
provides a criterion for the integrity of the method. The performance 
of each method in a regional, as well as a global sense is examined, and 
finally the results of the three methods are intercompared, taking into 
account also the c:>mplexity (or simplicity) of each method, its effi- 
ciency and the computational time that it requires. 
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Fig. 4 


Range versus time graph for the sparse data set from 
station 7943. 
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Fig. 5 Range versus time graph for the dense data set from 
station 7115. 
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3 . 2 . 3 . 1 Least squares approximation using monomials and 

Chebyshev polynomials as base functions . The existence and uniqueness 
of an n^^-degree polynomial that takes given values on a set of n + 1 
points of a closed interval [a, b] is guaranteed by the well-known 
theorem of polynomial interpolation [Davis, 1975]. The representation 
theorems of interpolation theory provide the tools for determining t'ns 
polynomial in various ways. Lagrange's and Newton's formulae are the 
more often quoted solutions for this problem. It is equally well known, 
though, that as the number of given points increases (the degree of the 
polynomial increases, too) insurmountable problems arise from the compu- 
tational aspect of this polynomial. Furthermore, even if the numerical 
problems could be overcome, the resulting polynomial will exhibit such 
strong oscillations between the fiducial points that it would be imp'>s- 
sible to use it as a reliable approximation in those intervals. Despite, 
therefore, Weierstrass uniform approximation theorem and Walsh simul- 
taneous interpolation and uniform approximation theorem, in practice we 
must find a working alternative free of the aforementioned drawbacks. 

One such alternative that we study in this section is the appl. cation of 
"best approximation," best in some sense soon to be defined. 

A natural requirement for any type of approximant to a function f 
is that it should be "close" to f. As soon as we define the notion of 
"closeness" quantitatively, we have established a criterion for deter- 
mining the best approximant of f in that sense. For the specific prob- 
lem encountered in this investigation, we can restrict ourselves to the 
theory as applied for normed linear spaces of finite dimension [Davis, 
1975]. We chose to work with the 2-norm, since in that case its 
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interpretation as the geometric distance between two elements of the 
space is simple and intuitively appealing. For an element x of the 
n-dimensional real space the 2-norm is defined as 



where x.j (i=l, .... n) are the components of x. If x,y are two elements 
of a metric space R*^, and d denotes the distance function in R*^, then 
enforcing 

d(x, y) = II x-y II (34) 

to hold for all x,yER^, makes R*^ a normed linear space. 

An alternate and more appropriate way of obtaining a normed linear 
space is to start with an inner product space. The inner product is a 
"two-slot machine" similar to the distance function, possessing linear- 
ity, symmetry, homogeneity and positivity, and denoted by <•, •>. If 
we force in an inner product space the following equality to hold 

II X II = /<x, x> (35) 

then we obtain a normed linear space. Depending on the definition of 
the inner product, we obtain different norms for the resulting normed 
linear space. As they are all results of the more primitive concept of 
an inner product, we call them "induced norms." 

We have already chosen to work with the 2-norm, and it is easy to 
see from (33) and (35) that the inner product should be defined as 

n 

<x, x> = X] (36) 

i = l ^ 
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to be consistent with the rest of our formulation. We can now define 
the best approximation in terms of closeness under the 2-norm: 

Definition : The best approximation of y by a linear combi- 

nation of (xi, Xj^) is aiXi + ... + aj^Xj^ if the following 
inequality holds for every choice of the constants Ai, ..., Aj^: 
e = II y - (aiXi + ... + a|^X|^)|| ^ || y - (AiXi + ... + Aj^X|^)|| (37) 

The quantity e is called the error norm, and it is obvious from the 

above definition that the best approximation is the one which minimizes 

e. Since the Xi, .... span a subspacc R , which contaiiis the approx- 

imant a^Xi + ... + aj^Xj^, the above inequality can be illustrated as a 

1 / 

projection of y onto this subspace R . The error e then can be viev/ed 
as the perpendicular distance from y to R , and the approximant as the 
projection of y onto R , as shown in Fig. 6. 



Fig. 6 Geometric interpretation of the best linear approximation. 

It can be shown that for a given inner product, solving the pro- 
jection problem is equivalent to solving the best approximation problem 
If the solution to either problem exists, then it is unique. Further- 
more, we can guarantee the existence of the solution (therefore its 
uniqueness, too) by choosing x,, ..., Xj^ to be linearly independent and 
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k >. n. We can now construct the projection of y onto R and thereby 
obtain at the same time the best approximation under the 2-norm, better 
known as least squares approximation. 

Let Xi, X|^ be a spanning set for R . By means of the Gram- 

Schmidt orthogonal ization process we can form an orthogonal basis ei, 

k k 

for R . We are seeking the projection s of y onto R so that 

k k 

(y-s) _L R or equivalently (y-s)_[e. for all e^. in R . Since {e^} are a 

k k 

basis for R and s is an element of R , we can find {c^. } constants so 
that 

k 

s = Z c. e. (37) 

i=l ^ ^ 


The error vector y-s being orthogonal to all e^. satisfies 


or 


< y-s, e^ > = 0 i = l, . . . , k 

< y. e. > = < s, e^. > i = l, . . . , k 


and by (37): 

< y, e. > = c. 11 e.p 

which leads to 

_ < y. e. > 



(38) 

(39) 


(40) 


(41) 


The projection s is therefore fully determined and by the equivalence 
theorem the representation (37) is the best approximation of y in the 
least squa*^es sense. 

If we v;ant to determine s in terms of the original spanning set 
{x^. }, then we modify (38), (39) and (40) accordingly: 
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1 


f -» 


■ • g, ' . , \ 

< y-s, > = 0 1 = 1 , .... k 

< y, X. > = < s, x^. > i=l, . . . , k 


and since now 


we get 


k 

S = E c, X, 

k 

< y. X, > = E c < X , X, > 

' X=1 


1 * 1 . 


(42) 

(43) 


(44) 

(45) 


The last equation can be written in matrix form as 

,^[<y, x.>]i = [Cj]^ 

where 


(46) 


(47) 


and it can be easily verified that the matrix G is symmetric. Equation 
(46) is the normal equation of the problem. It is also known as the 
Grammian of {x^. } with respect to the adopted inner product. The solu- 
tion of (46) yie''ds the sought for coefficients {c^l and the required 
inverse of G is guaranteed by the linear independence of the basis 
elements {x^. }. 

Two different sets of basis elements (x^l have been used in this 
study. The set of monomials 

M(t) = t^ )l=0, ..., k-1 (48) 

and the Chebyshev polynomials 4„(t) defined as 


To(t) = 1 
Ti(t) = t 
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T^(t) = 2t T^_^(t) - T^_2(t) i=2. .... k-1 (49) 

The former result in the widely used "polynomial fit" to the given data, 
while the latter, due to certain properties they have, show a superior- 
ity from the computational point of view that becomes increasingly more 
evident as the degree of approximation (the dimension of R , k) in- 
creases. In theory the two fits (for the same degree k) should be equi- 
valent, since T(t) are linear combinations of M(t) and vice versa, and 
therefore span the same space R . The former have the interesting prop- 
erty of being orthogonal with respect to summation though, when the 
summation is carried over specific points in the interval [-1, 1]: 

tn f 0 . i f 3 

[ m+1 , i = j = 0 

where {tj^} are the zeros of defined as 

^ (11+1)2 

This property is the basis for the excellent from all aspects Chebyshev 
interpolation [Dahlquist and Bjbrk, 1974; Snyder, 1966], when we can 
choose the distribution of the given data. It seems, however, that for 
a dense distribution of data sums of products of T^(t)'s, such as those 
encountered in the formation of the normal equations (46), tend to be- 
have in much the same way. The growth of the off-diagonal elements of 
the G matrix is slower in this case (with respect to the degree k) com- 
pared to a fit with monomials Mj(t). Furthermore, the loss of signifi- 
cance (due to the finite nature of the computer) is much less serious 
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here than it is for the classical monomial fit. We can thus use higher- 
degree Chebyshev expansions than polynomial ones without the risk of 
divergence and thereby absurd results due to an ill-conditioned normal 
matrix. 

The question then arises as to what should be the choice of the 
degree k. Obviously, for a convergent fit the error will decrease as 
the degree increases. The oscillation of the resulting approximant 
will also increase though, and we must find a way of stopping at some 
optimum degree before that and the loss of significance make the approx- 
imant worthless. Since some of the criteria that we have studied are 
based on statistical concepts, we first have to introduce such concepts 
to our approximation process, which so far is of a purely mathematical 
nature. The simplest way of doing this is to modify our definition of 
the inner product (36) by including a weight function w(x^) which is 
assumed to be positive definite: 

n 

< X, X > = H w(x.) X.2 (52) 

w 1 1 

The meaning of w(x^) can be that of assigning various degrees of 
importance to each x^ in a relative sense. One natural choice in the 
least squares approximation is to gauge importance against the amount of 
information contained in each observation for the parameters of the 
problem. To do this we must have some measure of the level of observa- 
tional errors, and this can be achieved statistically by defining their 
distribution. If the distribution is the normal, which is usually the 
case, then it is fully defined through its first two moments, the mean 
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u and the variance o^. This being the case, it can be shown [Rao, 1973] 
that the Fisher information measure of each observation on the y param- 
eter of its distribution is 



(53) 


J’ iu) is a function that complies with the requirements of the weighting 
function, and from its definition it seems that it is one fit for this 
purpose. The weights are therefore determined as 


w(x.) = i (54) 

' o.^ 

1 


where in our case we have used o^. ^ (i = l, •••» n), since we have no 
means to discriminate between the observations. In matrix notation, 
the inner product can now be expressed by the following quadratic form; 

<x,x>=x^Wx (55) 

w 

where 

W = diag [w(x^)] (56) 


If one can further make the assumption that the errors are dis- 
tributed independently, then W is the inverse of their variance matrix 
L, and in this case the approximation is more appropriately called 
"minimum variance estimation" [Rao, 1973]. We can now, in light of the 
above discussion, speak of statistical properties of the approxiniant, 
and in fact we should also change our terminology, substituting estima- 
tion for approximation, estimator for approximant, and residuals for 
errors. This generalization of the process of approximation opens the 
way for use of various tests of significance, devised and applied in the 
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theory of linear statistical inference [Rao, 1973; Pollard, 1977]. In 
the tests that we have performed with the two selected sets of data we 
have looked at possibilities such as the convergence of the root mean 
square error of the fit (rms of fit), the rms of the recovered values 
at the ground truth points, and the significance of the change of the 
weighted sum of squares of the residuals between successive fits. 

3. 2. 3. 2 Interpolation with cubic spline functions . Despite 
the fact that theoretically one can always determine an interpolatory 
polynomial for the station-satellite range function, we have already 
seen that the large number of data points makes its computational aspects 
awkward and its qualitative and quantitative value questionable. To 
avoid the erratic behavior of the polynomial in between data points, we 
must keep its degree low. If, however, low-degree polynomials are 
fitted to the data, then we must tolerate the fact that they only approx- 
imate the function, that they do not reproduce the function at the fidu- 
cial points, and additionally, that some filtering of the high frequen- 
cies in the data is unavoidable. An alternative that exhibits the 
described low-degree-polynomial behavior and tne reproducibility of the 
function at the fiducial points is the use of spline functions [Spath, 
1974]. 

Spline functions (SF's) are curves consisting of low-degree poly- 
nomials each of which is defined over the closed interval [t^, t^.^^] of 
two successive fiducial points. These elements of the SF are connected 
at these nodal points in such a way that their derivatives (to maximal 
order) exist. Since the SF passes through all nodal points, it 
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reproduces the interpolated function at these points exactly; and fur- 
thermore the use of low-degree polynomials suppresses the undesired 
oscillations between these nodal points. Provided that some boundary 
conditions are specified for the behavior of the SF at the endpoints of 
the interpolated interval, the existence and uniqueness of the SF is 
guaranteed [Spath, 1974]. 

As is the case with least squares approximation, the boundary con- 
ditions for the SF determine an "optimality" criterion which the result- 
ing SF satisfies. A natural cubic spline s(t), tG[a,b], for instance, 
is determined with the boundary condition 

s"(a) = s"(b) = 0 (57) 

and it is shown in [Spath, 1974] that s(t) minimizes the following 
integral 

b 

1(f) = I [f"(t)]^dt (58) 

a 

where f(t) is the interpolated function. Equivalently, one can say that 
the above SF is the solution to the variational problem 

b 

I(f) - / [f"(t))]^dt = a minimum (59) 

a 

with the aforementioned ' jundary conditions supplemented by the addi- 
tional constraints that the resulting function passes through the given 
nodal points. Even though splines are interpolatory in nature, with a 
reformulation of the problem they can be "fitted" to the data in some 
optimal sense--most naturally the least squares sense. So one can con- 
ceivably solve the problem of approximation and filtering simultaneously 


63 



ORICK-JAL PAGir f3 
OF POOR QUALirV 


with this option [Spath, 1974]. Since in our case the weighting func- 
tion for the observed nodal values is not known a priori, we refrain 
from using smoothing splines. In what follows, we summarize the compu- 
tational aspects of cubic splines as we have applied them in this inves- 
tigation. 

We are given the values of the station-satellite range r^ at n 
epochs, ti<t 2 <. . .<t^. We seek to determine the constituents f^. (i = l, 
.... n-1) of the cubic spline s(t), each f^ being a cubic polynomial, 
under the constraint that s(t) is twice differentiable at the n epochs 
ti, ..., t^. Each of the f^'s is defined on the corresponding range 
[ti* t.^j]. We adopt the following form for f ^ ; 

f,(t) = (t-t,)= a * a,_3 (t-t.) + a,^^ (60) 

so we will have four unknown coefficients for each of the n-1 f/s, 
4(n-l) total. To exploit the constraints for derivative continuity and 
existence, let us first establish the following notation for brevity: 


at, 

= ‘ui ■ S' 

Ar . 
1 

•"i+l " 

r . 
1 

= f,(t,) 

•"i+l 


r: 

1 

II 

•^i+l 


r'.' 

1 

= f. (t.) 
1 ' 1 

•"i+l 



We can now write the following set of equations for each of the n-1 
intervals [t., t.^^]: 


r . = a . - 

1 1 ,4 


•■hi = * ^.2 "‘i ^ ^,3 “i * ^,4 


(62) 

(63) 
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r! = a. 


r ' 

*^1+1 


i ,3 

3a. , At? + 2a. o At. + a. 
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i,l “"i 


1,2 1 1,3 


= 2a. 


i+1 


i,2 

6a . , At. + 2a 
1 » 1 1 


i,2 


(64) 

(65) 

( 66 ) 

(67) 


After some algebra, (62), (63), (66), and (67) can be solved for 


the four unknown coefficients (a.:.:. J=l> •••» 4}: 

* J 


^.1 = wt: '•'m - '■P 


^,2 " 2 "'i 


a y, = r. 
1 ,4 1 


( 68 ) 

(69) 

(70) 

(71) 


Substituting (68) - (71) into equations (64) and (65) we obtain 
for the first derivatives 


and 


-■i 


Ar , , 

r' = -r^— + r At , (2r" + r" ,) 
n At 1 6 n-1 n n-1' 

n-1 


(72) 


(73) 


Fran the continuity constraints for the first derivative we must 
force the following equality to hold at the nodes; 


n <‘itl> = ^ial “itl) 
Substituting f-^om (64) and (65) we get 
of the coefficients (aat): 

* sJ 


i=l, . . . , n-2 (74) 

the explicit constraint in terms 
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( 75 ) 


We now substitute from (68) - (71) for the (a-J in terms of the second 

* V 

derivatives rV, etc. and obtain the following working expression: 


4t,r;' t (2st. + = 6 




At, 


i = 1, 


i+1 

n-2 


Ar . 



At, 


(76) 


Obviously this expression cannot be written fcr the first and last 
nodal points, for there are neither prior (for the first) nor posterior 
(for the last) information on which the estimation of the r'{ and r^ can 
be based. We therefore need to provide this information through the 
boundary conditions. In this case the first and last equations of the 
(76) system are modified appropriately. When written in matrix form, 
(76) is as follows: 


2(Ati+At2) 

A t2 

0 

>2’ 

A t-2 

2(At2'*’At3) At3 


r'3 


At 3 

^^n-2 

r" 

n-2 

0 

Osi 

1 

c 

<3 


r" , 

L n-lj 




Ar, 


Ar, 


Ata 

Atti ■ ' 

Arj 


A 1 3 

Atz / 

Ar o 

Ar T 

n-2 

n-3 

At o 

’ At 0 

n-2 

n-3 


Ar „ 

n-2 

"‘n-1 

^^^n-2 


At 


n-1 n 


( 77 ) 
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For the system (77) to have a unique solution the coefficient 
matrix must be invertible. This matrix is symmetric, tridiagonal, with 
positive diagonal elements, and diagonally dominant. None of the pri- 
mary determinants vanish therefore, and in fact they are all positive. 
The matrix is thus positive definite, and its Cayley inverse exists 
(unique inverse). For the particular case of SLR data interpolation, 
the second derivatives at the two end points are not known explicitly, 
and so we must look for an alternative set of boundary conditions. 

Since the observed range is a slowly varying function, we can safely 
assume that r'i and r^ are related linearly to ra and r^ ^ respectively. 


r'l = 0 rj' 


r" = T r", 
n n-1 


(78) 


where o, t are constants which we must choose. The new set of boundary 
conditions (78) requires that the first and last of tlie equations (77) 
be changed to 


(2+o)Ati + 2Ata Ata 


Atn_2 2At^_2 + (2+i)At|^_j 


/ (ill Ari \ 

■ 4t, ) 


ri 


n-1 


(79) 
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Finally, the choice we have made for o and t is to set them 
equal to unity. The rationale behind this is that the derivatives of a 
polynomial become in general smoother as their order increases, and since 
we have a smooth function to begin with, there cannot possibly be a sig- 
nificant change in its second derivative over the few seconds that the 
two nodes are apart. Besides this, as shown in [Spath, 1974], the vari- 
ation of the boundary conditions results in insignificant variations in 
the interior intervals and therefore a stable interpolation spline will 
be obtained for any reasonable choice of these conditions. The final 
form of (79) is therefore (o=x=l): 



3Ati+2At2 
A t2 


0 


A t2 

2(At2+At3) 


At 3 
** • 


At 


n-2 


1 

o 


'r-; ■ 




• 


• 



• 



r" . 
n-1 


Ar2 

At2 


Ari 
At 1 


Ar; 

At; 


Ar2 

At2 


Ar 

At 


n-1 

n-1 


^'"n-2 


At 


n-2 


(80) 


Solution of (80) yields the required second derivatives rV (i=2, 

..., n-1) for the determination of the 4(n-l) coefficients (a-.) 

^ J 

through (68) - (71). At this point the spline s(t) is fully determined 
by the n-1 cubic polynomials f^. (t) (i-1, ..., n-1) with the general ex- 
pression (60). Values of the interpolated function can be be obtained 
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the ground truth points to get an 
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We can use the interpolated values at 
idea of how well s(t) approximates 


3 , 2 . 3 . 3 Comparison of least squares estimates and cubic spline 
interpolants for the range ^unction . We present and compare in this 
section the results obtained from a number of tests we performed with 
the SLR data collected during two passes of satellite Lageos over the 
Australian station at Orroral (7943) and the U.S. Station at Goldstone 
(7115). In all cases, except foi^ spline interpolation, the domain of 
the function has been transformed to [-1, 1] through the following 
equation 


T(t) = 


2t - (tg + tg) 


(81) 


where U, tr are the actual epochs of the first and last observations 

b L 

available. The above transformation greatly improves the conditioning 
of the normal equations, since raising t to high powers can now result 
only in losing some significance if t is very small. If t could be 
larger than one, then we ran the risk of exceeding the exponent magni- 
tude limit of our computer. 

With reference to the conditioning of the normals, we have also 
found helpful the scaling of the matrix in what is commonly termed 
"correlation form" in statistics. This technique results in making 
uniform the range of values of the elements of the matrix, thereby im- 
proving its numerical properties. The scaling is done jy means of the 
square roots of the diagonal elements. We denote the original system 
of normal equations by 
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Nx = U 


(82) 


Then let 

C = diag (83) 

with |C| ^ 0 since N is positive definite. We put 

N = CNC (84) 

and therefore 

N"^ = C’^N'^C'^ (85) 

Pre- and post-multiplying (85) by C we obtain 

CN‘^C = CC‘^N'^C"^C (86) 

or 

N"^ = CN’^C (87) 

from which we obtain in combination with (82) 

X = CN"^C U (88) 

Although this procedure does not alleviate numerical problems com- 
pletely, it seems to improve the solution, especially in the case of a 
good distribution of data, as can be verified from Table 3. When the 
base functions are the Chebyshev polynomials, then the original normals 
are already uniform, and we see no change in the results of the two 
sol utions. 

The quality of the approximation depends not only on the global 
behavior of the approximant, but on the local as well. A generally good 
fit therefore can give poor results in some regions where the data dis- 
tribution is worr? than in the rest of the data set. When the differ- 
ences oetween the o^'served ranges and the approximated at the ground 
truth points are examined, it is helpful to know whether we are working 
in a dense region or a sparse one. We have determined the time intervals 
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Table 3 Comparison of Monomial Fits Using N and N 


Tes 1 

Case 


Dens 

e U a t 

a 

s 

p ii r ^ 

e Dal 

a 

Ue (tree 

Nurimi 1 

HflS 

Me a 11 

lie c o ve I* y 

cru 

IIMS 

Me a n 

lie CO ve ry 

Cl’U 

o f 

Kcfiis 

FI t 

Error 

IIMS 

T 1 me 

FI t 

Erro r 

1U'L> 

T 1 me 

Fit 

Used 

( m) 

( ni) 

< iti) 

(s) 

( III) 

(m) 

( in) 

( s ) 


N 

0.24 

0.01 

0. 12 

1.3 

3.5 

-0.60 

2.36 

0. 13 

IG 











N 

0.22 

-0.01 

0.00 

1.3 

0.5 

0. 10 

0.59 

0. 15 


N 

7.21 

0.23 

6.37 

1.7 

16 19.9 

-406.20 

1570.30 

0. 10 

19 











N 

7.39 

-0,50 

4.61 

1.7 

1490.9 

0.25 

14.04 

0. 10 


ORICrJ.^l PAGE T5 
OF POOR QUALiTY 


Table 4 Distribution of Ground Truth Points, 
Sparse Data 


Station No. : 7943 
Observation Epoch 

< 

T 

Tl X — 

> 

YYMMDD 

HHimSS.SSSS 



T 

800814 

1 10626.8209 

7.50 

15.00 

22.50 

8008)4 

1 10719.3610 

37.54 

22.53 

60.07 

800814 

110842.0612 

7.50 

15.00 

22.50 

800814 

1 10934.5613 

15.00 

7.50 

22.50 

800814 

111019.5615 

15.02 

7,50 

22.52 

809814 

nil 4.5216 

22.46 

22.46 

44.92 

800814 

1 1 1 134.4617 

7.48 

7.49 

14.97 

800814 

1 1 12 4.4618 

15.00 

7.50 

22.50 

800814 

1 1 1219.4319 

7,47 

15.00 

22 47 

800814 

1 1 1256.9021 

7.47 

7.50 

14.97 

800814 

1 1 1319.4022 

15.00 

7.50 

22.50 

800814 

1 1 1356.8823 

22.48 

15.00 

37.48 

800814 

1 1 1526.8828 

15.00 

45.02 

60.02 

800814 

1 1 1626.9031 

7.50 

59.96 

67.46 

800814 

1 12019. 1846 

7,50 

15.00 

22.50 


7i 


to the nearest observation prior to and after each of the ground tru.h 
points for both data sets. The epoch of these observations and these 
time intervals are displayed in Tables 4 and 5, where T^ is the interval 
in seconds to the nearest preceding observation, and T^ is the interval 
to the nearest succeeding one. 

In all tests we have used expansions starting from degree five all 
the way up to 20. The summaries which are discussed here display only 
partial results which are enough to indicate the performance of each 
method. From Table 6 it is obvious that even for the dense pass the 
approximation with monomials starts diverging after degree 16. This is 
also true for the other pass, as ca" 'een from Table 7. The results 
for solutions up to degree 14 are . for monomials and Chebyshev 

polynomials. From thereon though, wf.i.e the former diverge, the latter 
converge with no major problems. The results given in Tables 8 and 9 
Ccn be compared to those of Tabels 6 and 7 respectively to verify this. 
Only the recovered range data quality becomes poorer at higher degrees 
for the sparse data set, and this is caused by the distribution of the 
data, rather than instabilities in approximation technique. The denser 
pass shows a very good recovery even at those high degrees. 

The computational time is about ten percent higher for the Cheby- 
shev solutions compared to the monomials, but the increased stability 
and quality of the solution seems to be well worth it. It is therefore 
recommended that if one has to choose between these two solutions one 
should always go with the Chebyshev expansion rather than the monomials. 
The question that arises next is how to decide which is the lowest degree 
that gives a satisfactory representation of the data. Using the 
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Table 5 

Distribution 
Dense Data 

of Ground 

Truth 

Points, 

Sta t Ion No . : 7113 

< 

— T - 

> 

Observation Epoch 

< Tl X— T|^”> 

YYMMDD 

HHfOtSS.SSSS 


'Tr 

T 

000814 

131626.0577 

1.00 

1.00 

2.00 

BOOB 14 

131653.0577 

1.00 

1.00 

2.00 

800814 

131714.0377 

1.00 

1.00 

2.00 

800814 

131735.0577 

1.00 

1.00 

2.00 

BOOB 14 

1318 2.0577 

1.00 

1.00 

2.00 

BOOB 14 

1320 9.0577 

1.00 

1.00 

2.00 

800814 

132032.0578 

1.00 

1.00 

2.00 

800814 

132052.0578 

1.00 

1.00 

2.00 

800814 

132113.0578 

1.00 

1.00 

2.00 

800814 

132134.0579 

1.00 

1.00 

2.00 

800814 

132156.0579 

1.00 

1.00 

2.00 

8008 14 

132217.0580 

1.00 

2.00 

3.00 

800814 

132240.0581 

1.00 

2.00 

3.00 

800814 

1323 4.0581 

1.00 

1.00 

2.00 

800814 

132325.0382 

1.00 

1.00 

2.00 

800814 

132345.0583 

1.00 

1.00 

2.00 

800814 

1324 9.0584 

2.00 

1.00 

3.00 

BOOB 14 

132436.0585 

3.00 

2.00 

5.00 

800814 

132518.0586 

3.00 

1.00 

4.00 

800814 

132347.0588 

3.00 

1 .00 

4.00 

B 00814 

132623.0589 

1 .00 

1.00 

2.00 

800814 

132652.0591 

1.00 

1.00 

2.00 

860814 

132720.0392 

6.00 

1.00 

7.00 

800814 

132745.0394 

2.00 

1.00 

3.00 

800814 

13281 1.0595 

2.00 

1.00 

3.00 

800814 

132833.0597 

1.00 

1.00 

2.00 

800814 

132854. 0508 

1.00 

1 .00 

2.00 

8006 14 

13291 , 5.0599 

1.00 

1.00 

2.00 

800814 

132039.0601 

1.00 

I .00 

2.00 

800814 

1330 4.0602 

1.00 

2.00 

3.00 

800814 

133038.0604 

1.00 

1.00 

2.00 

800814 

1331 5.0606 

1.00 

1.00 

2.00 

BOOB 14 

133137.0609 

3.00 

1.00 

4.00 

800814 

1332 9.0611 

2.00 

1.00 

3.00 

800814 

133247.0614 

7.00 

3.00 

10.00 

800814 

133329.0617 

2.00 

2.00 

4.00 

800814 

133427.0622 

4.00 

1.00 

5.00 

800814 

133528.0627 

14.00 

14.00 

28.00 
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Table 6 Least Squares Approximation with Monomials, Dense Data 


Station: 7115 Date: 800814 Obs. Approx.: 754 Obs. Recov.: 38 


— 

Degree 
of Fit 


Average 
Recovery 
Error (m) 

RMS of 
Recovery 
Errors (m) 

Average Serial 
Correlation for 
Recovered Obs. 

1 

CPU 

Time (ms) 

6 

5.373 

-0.477 

5.316 

0.947 

353 

1 

8 

0.249 

-0.014 

0.130 

0.923 

496 ! 

10 

0.220 

-0.016 

0.079 

0.893 

634 

12 

0.220 

-0.014 

0.084 

0.860 

803 I 

14 

0.219 

-0.013 

0.089 

0.818 

1064 

15 

0.219 

-0.015 

0.081 

0.801 

1177 1 

16 

0.220 


0.084 

0.785 

1315 1 

1 

17 

0.323 

0.024 

0.133 

0.770 

1436 ' 

j 

18 

0.505 

-0.162 

0.334 

0.751 

1581 

20 

35.919 

-1.033 

9.041 



0.709 

1859 



Table 7 

Least Squares Approximation with Monomials, Sparse Data 

Station 

: 7943 Date: 800814 

Obs. Approx 

. : 41 Obs. Recov 

. : 15 i 

1 

Degree 
of Fit 


Average 
Recovery 
Error (m) 

RMS of 
Recovery 
Errors (m) 
— 

Average Serial 
Correlation for 
Recovered Obs. 


6 

1.350 

0.110 


0.859 

43 

8 

0.381 

0.082 

0.458 

0.806 

57 

10 

0.340 

0.051 

0.472 

0.7^7 

70 

12 

0.339 

0.045 

0.473 

0.730 

91 

14 

0.336 

0.032 

0.474 

0.645 

121 


0.332 

_.082 

0.561 

0.573 

135 i 

16 

0.460 

0.096 

0.592 

0.484 

mm 

17 

18.399 

- 1.435 

5.586 

0.412 

166 j 

18 

272.217 

-57.402 

222.250 

0.355 

182 

1 

ro 

O 

68626.4 

>10^ 

48583.9 

0.277 

210 

1 
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Table 8 Least Squares Approximation with Chebyshev Polynomials, 
Dense Data 


Station 

: 7115 Date: 800814 

Obs. Approx 

. : 754 Obs. Recov.: 38 j 

Degree 
of Fit 

RMS of 
Fit (m) 



Average 
Recovery 
Error (m) 

RMS of 
Recovery 
Errors (m) 

Average Serial 
Correlation for 
Recovered Obs. 

CPU 

Time (ms) 

mm 

mmm 

-0.016 


0.801 

1296 

H 

RMI 

-0.014 

0.086 

0.785 

1433 1 

17 

0.219 

-0.015 

0.084 

0.770 

1580 i 

1 

18 

0.218 

-0.015 


0.751 

1731 1 

1 

19 

0.218 

-0.015 


0.730 

1898 

20 

0.218 

-0.015 

0.087 

0.709 

2064 

• 

- -- 


Table 9 Least Squares Approximation with Chebyshev Polynomials, 
Sparse Data 


Stdti on 

7943 Dace: 800814 

Ols. Approx. 

41 Obs. Recov. 

: 15 ; 

Degree 
of Fit 

RMS of 
Fit (m) 

Average 
Recovery 
Error (m) 

RMS of 
Recovery 
Errors (m) 

t 

Average Serial 
Correlation for 

Recovered Obs. 

— 

1 

CPU ; 
Time (ms)i 

I 

1 

15 

^;.3 

0.070 

0.531 

.... 

0.573 

1 

140 

16 

0.319 

0.015 

0.458 

0.484 

156 1 

17 

0.319 

0.011 

0.455 

0.412 

173 1 

18 

0.318 

0.386 

1.592 

0.355 

191 : 

19 

0.318 

0.642 

2.561 

0.294 

209 1 


0.317 

-0.762 



2.964 

0.276 

230 

1 


Note; CPU times above refer to the Amdahl V8 computer of The Ohio 
State University Instruction and Research Computer Center. 
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statistics on the recovered ground truth data gives some indication of 
the quality of the fit. It is, however, conditional on the selection, 
distribution and number of such points within the available data. It is, 
therefore, a very local test and does not provide a measure of the glob- 
al performance of the estimator. The weighted sum of squares of the 
residuals, on the other hand, has this property, and it is a standard 
procedure in this type of problem to test the significance of the change 
of this statistic between solutions based on the same data. Hamilton 
[1964] derives the R-test to test the null hypothesis 

Hq: the degree fit is as good as the (k+1)^^ degree fit 
based on the variance ratio F-test. 

The test statistic is the following: 

(V^PV)^, 

R = (89) 

(V^PV)^,1 


where the subscript indicates the fit from which the residual norm has 
been computed. The theoreticdl value is obtained as 


^^^1, (n-k-1), a ■ n-k-1 ^1, n-k-1, a 


+ 1 


(90) 


with significance level 100a%. The hypothesis is rejected at this level 
if R The test must be used carefully, and one should make certain 

that the process has reached a stable convergence before one applies the 
test. Furthermore, the hypothesis to be tested should actually be a 
"chain-type" one, in the sense that we test from the highest degree down 
to the lowerst one where we find that we must reject the hypothesis. 

For example, if we test in sequence 
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Hoi: a 19th-degree fit is as good as a 20th 

H02: an 18th-degree fit is as good as a 19th 

H07: 3 13th-degree fit is as good as a 14th 

and we find that Ho? is rejected while all previous hypotheses were 
supported by our results, then we can say that a 14th-degree fit is as 
good as any of the higher-degree fits at the tested significance level 
Carrying out this test for the two test data sets, we have found that 
for the dense data set the lowest degree acceptable fit is for degree 
ten, and for the sparse data set it is degree seven. The test statis- 
tics are given in Table 10. 


Table 10 R-Ratio Test Results (a = 0.01) 


Data 

Set 

— 

Dense 

Data 


1 

Sparse Data 

Degree 

T 

v'pv 


V[,99% 

Ho 

V^PV 



Ho 

5 

232.05 

1.0661 

1.0080 

R 

2.6828 

3.5881 

1.329 

R 

6 

217.67 

72.036 


R 

0.7477 

11.2946 


R 

7 

3.0217 

6.4580 


R 

0.0662 

1.1107 


1 

A 

8 

0.4679 

1.2494 


R 

0.0596 

1.0700 


A 

9 

0.3745 

1.0272 


R 

0.0557 

1.1740 

1.252 

A 

10 

0.3646 

1.0006 


A 

0.0474 

1.0023 


A 

11 

0.3644 

1.0023 

1.0080 

A 

0.0473 

1.0040 

1.230 

A 

12 

0.3635 




0.0471 





Note: A - Accept, R - Reject 
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We have finally examined the statistical independence of the esti- 
mated range values at the ground truth points. Even if we assume that 
the available data are contaminated with uncorrelated random errors, 
the recovered ranges have approximation errors that are quite strongly 
correlated. The closer these points are to each other, the higher the 
correlation between every successive pair. The level of correlation 
drops as the degree of the fit is increased or as the data set becomes 
sparser. We have computed the average value of this serial correlation 
for each of the fits, and it is listed along with the other statistics 
in Tables 6 and 7 for fits with monomials, and Tables 8 and 9 for fits 
with Chebyshev polynomials. The choice of base functions has no effect 
whatsoever on this correlation. For the two fits that have been selec- 
ted on the basis of the test, the corresponding correlation levels 
are 83% for the seventh-degree fit on the sparse data set and 89% for 
the lOth-degree fit on the dense data set. Since a significant level 
of correlation is seen even between every fourth or fifth observation, 
one would have to include a full weight matrix in any subsequent use of 
these ranges if a meaningful result is sought. This, however, would 
make the use of such data very cumbersome and increase the computational 
effort beyond reason. 

A way to circumvent this difficulty, without compromising on the 
assumed statistics of the estimates, is to use the cubic spline interpo- 
lation. Results for the two data sets are given in Tables 11 and 12. 
Since the splines fit exactly at the data points, the rms of the fit is 
identically zero and cannot be used as an indicator of the quality of 
the fit. We can use though the recovered ground truth data statistics 
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Table 11 Interpolation with Cubic Splines - Dense Data Sets 


Test 

Observ. 

Approx. 

Observa- 

tions 

Recovered 

Average 
Recovery 
Error (m) 

RMS of 
Recovery 
Errors (m) 

CPU 

Time (ms) 

I 

754 

38 

-0.032 

0.210 

125 

II 

754 

37 

-0.007 

0.147 

125 

Table 12 Interpolation with Cubic Splines 

- Sparse Data 

Set 

n 

Test 


Observa- 

tions 

Recovered 

1 

Average 
Recovery 
Error (m) 

RMS of 
Recovery 
Errors (m) 

CPU 

Time (ms) 

I 

41 

15 

-0.202 

0.944 

25 


to compare with the least squares estimates. Because splines are very 
sensitive to data distribution, we have recomputed the statistics for 
the dense pass after deleting the last entry, which as can be seen from 
Table 5 is isolated and does not conform with the rest of the test 
points. This results in a significant improvement of the statistics. 

One should justify this in the sense that uniformity must exist if one 
wants to obtain results of seme quality. At any rate, it is apparent 
that the results here show a higher rms error by a factor of two com- 
pared to the least squares estimates, but at the same time the computa- 
tional effort is about three times smaller. In addition to that, be- 
cause the interpolation over each interval is done using a different 
constituent of the spline, interpolates from different intervals are 
totally uncorrelated. One way to show this is the following experiment. 
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We have recomputed the spline curve for the dense data set, only this 
time we truncated the data set to data collected prior to 13^32^00^, 
and thus we are not recovering the last five ranges of the ground truth 
data. If this is done for the least snuares estimator, the errors for 
the recovered ranges change, because the fit is done using information 
from all data points simultaneously. In the case of the spline curve 
though, only local information is used, and as can be seen from Table 
13 the errors of recovery are identical at the common points for both 
the complete (po - Sc(t)) and the shortened (po- Ss(t)) data sets. 

It can now be safely assumed that the interpolation errors for the 
cubic spline are uncorrelated, and *:herefore the error character! sties 
of the interpolated ranges are not altered by this p<^ocess. For a uni- 
form dense distribution of the base points in the region where the 
ranges are interpolated, the mean value of the recovery errors is below 
the centimeter level, and their fluctuation does not indicate a disper- 
sion significantly different from that of the original data (cf. Table 
11). It is very simple to check the quality of the interpolated ranges 
by examining the time 'ntervals T|^ and Tp. When tht;e intervals are 
significantly different from their average value over the entire pass, 
then it might be wise to delete that range from the data set or at 
least give i. a smaller weight in subsequent use. One cannot form 
strict ruler, to follow for this procedure; it is more easily solved on 
a case-by-case basis, and the action taken depends largely on the ex- 
perience and judgment of the investigator. For the purpose of this 
study, we feel that the right approach is to delete such ranges com- 
pletely. 
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Table 13 Comparison of Recovery Errors from Two Cubic Splines 
for the Complete and the Restricted Dense Data Sets 


Station No. : 7115 
Obaervatlon Epoch 

< 

Tl — -X — 

tr--> 

Recovery Error 

YYMMDD 

HHimss.ssss 

TL 

tr 

T 

( m> 

( m) 

BOOH 14 

131626.0577 

1.00 

1.00 

2.00 

0. 1000 

0 . 1 000 

BOOB 14 

131653.0577 

1.00 

1.00 

2.00 

0 1191 

0. 1191 

BOOB 14 

131714.0577 

1.00 

1.00 

2.00 

-0.0708 

-0.070B 

BOOB 14 

131735.0577 

1.00 

1.00 

2.00 

-0.0515 

-0.05 15 

BOOB 14 

13IB 2.0577 

1.00 

1.00 

2.00 

0. 1370 

0. 1370 

BOOB 14 

1320 9.0577 

1 .00 

1.00 

2.00 

O.OBIB 

O.OBIB 

BOOB 14 

132032. 057B 

1.00 

1.00 

2.00 

0 . 0B7B 

0 . 0H70 

BOOB 14 

1 32052. 057B 

1.00 

1.00 

2.00 

-0. 1059 

-0. 1059 

BOOH 14 

1321 I3.057B 

1.00 

1.00 

?.00 

-0. 1053 

-0. 10.53 

BOOB 14 

132134.0579 

1.00 

1.00 

2.00 

-0.0250 

-0.0250 

BOOB 14 

132156.0579 

1.00 

1.00 

2.00 

-0.0225 

-0.0225 

BOOB 14 

132217. 05B0 

1.00 

2.00 

3.00 

-0.0624 

-0.0624 

BOOH 14 

132240.0501 

1.00 

2.00 

3.00 

-0.0514 

-0.0514 

BOOB 14 

1323 4.05B1 

1.00 

1.00 

2.00 

-0.0791 

-0.0791 

BOOB 14 

1 32325. 05B2 

1.00 

1.00 

2.00 

-0. 12B9 

-0. 1209 

BOOB 14 

1 32345. 05B3 

1.00 

1.00 

2.00 

-0.0216 

-0.0216 

BOOB 14 

1324 9.05B4 

2.00 

1.00 

3.00 

-9. 1772 

-0. 1772 

BOOB 14 

1 32436. 05B5 

3.00 

2.00 

5.00 

. 0000 

0 . OOOB 

BOOB 14 

1325 IB.05B6 

3 . OO 


4.00 

e . 0304 

0 . U3U4 

BOOB 14 

1 32547. 05BB 

3.00 

1.00 

4.00 

0.0733 

0.0733 

BOOB 14 

132623. 05B9 

1.00 

1.00 

2.00 

0.2133 

0.2133 

BOOB 14 

132652.0591 

1.00 

1.00 

2.00 

0.0625 

0 . 0625 

BOOB 14 

132720.0592 

6.00 

1.00 

7.00 

-0.3168 

-0.3I6B 

BOOB 14 

132745.0594 

2.00 

1.00 

3.00 

-0. 1 194 

-O. 1 194 

BOOB 14 

I32B1 1.0595 

2.00 

1.00 

3.00 

-0.0176 

-0.0176 

BOOB 14 

132B33.0597 

1.00 

1.00 

2.00 

0. 1412 

0. 1412 

BOOB 14 

l32B54.ei;9B 

1.00 

1.00 

2.00 

-0.025B 

-0.0250 

BOOB 14 

1329 15.0599 

1.00 

1 .00 

2.00 

-0.2227 

-0 . 2227 

BOOB 1 4 

1.32939.0601 

1 .00 

1.00 

2.00 

0.0429 

0 . 0429 

BOOB 14 

1330 4.0602 

1.00 

2.00 

3.00 

0. 1007 

0. 1B07 

BOOB 14 

I3303B.0604 

1.00 

1.00 

2.00 

0. 1012 

0. 1012 

BOOB 1 4 

1331 5.0606 

1.00 

1 .00 

2.00 

-0. 1256 

-0. 1256 

HOOB14 

133137.0609 

3.00 

1.00 

4.00 

0.4690 

0.4690 

Bt.‘»ai4 

1332 9.0611 

2.00 

1.00 

3.00 

-0.0813 


O'^Oh ' 4 

133247.0614 

7.00 

3.00 

10.00 

-0 . 2007 

— 

tli>OG "‘t 

133329.0617 

2.00 

2.00 

4.00 

-0.0889 

— 

BOOB 14 

133427.0622 

4.00 

1.00 

5.00 

0 . 003 1 

— 

BOOB 14 

I3C52B.0627 

14.00 

14.00 

20. 00 

-0.9416 

— 



The ease by which spline interpolation can be applied in our 
problem and the properties that we have found it to posses encouraged 
us to use it as the most suitable method for interpolating the quasi- 
simultaneous ranges required for determining the station-satellite- 
station range differences. 
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4. THE ESTIMATION PROCESS 


4.1 Introduction 

The motions of the satellite and the observing stations in space 
have been described in the second chapter in terms of models that 
depend on a number of parameters. Some of these parameters are obtained 
from theories based on very long (time-wise) records of observations and 
therefore carry a great deal of confidence in them (e.g., precession, 
nutation). A considerable number of these parameters, though, are only 
approximately known, and part of the reason for requiring an adjustment 
of the observations is the improvement of their numerical values. The 
other more obvious reason, of course, is the smoothing of the errors 
inherent in every measurement process. In short, the adjustment process 
determines the unknown parameters based on the information contained in 
the discrepancies between the measured values of the observables and 
those computed from the assumed model. The operator that relates the 
corrections in the parameters to these discrepancies is the design 
matrix of the problem. In the usual case where there are redundant ob- 
serv Lions available, the row space of the design matrix has a dimension 
larger than that of its column space. Its rank then is determined by 
the dimension of its column space. If its columns are linearly indepen- 
dent, then the rank is equal to their number, the dimension of the 
column space, and the problem will have a unique solution. In the event. 
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though, that two or more of its columns are linearly dependent, the 
design matrix is rank deficient, its deficiency determined by the number 
of interdependent columns. When this happens, the problem does not 
admit a unique and unbiased estimate for the parameters. Special tech- 
niques must be employed in order to obtain even a unique estimate, 
preferably one with a minimum bias too [Rao, 1973]. Other techniques 
which are applicable in such a case, although with different properties 
from the previous one, have been reviewed in [Pavlis, 1979]. 

The fact that the relationship between observations and parameters 
may be a nonlinear one, as is the case here, further complicates the 
process. Although extensive literature exists for the linear model, 
that for the nonlinear one is very much restricted and hardly ever 
addresses the validity of extending statistical concepts established for 
the first for use in cases involving the second. Technically, the solu- 
tion is most uotally obtained by means of a Newton-Gauss iterative pro- 
cedure [Pope, 1974J. Starting from some approximation to the solution, 
the nonlinear relationship is expanded as a Taylor series retaining 
terms up to those including the first derivative. Assuming continuity 
and boundness for the original nonlinear function, the iteration might 
converge to the sought-for solution [Pope, 1972]. The convergence can 
be established by examining the percent change in the weighted sum of 
squares of the residuals between successive iterations. This test quan- 
tity can also be used to detect a divergent problem or an oscillating 
one. Possible explanations for these cases can be found, for instance, 
in [Hamilton, 1964; Pope, 1972; Uotila, 1975]. 
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Even when the problem has converged because of the nonlinearity, 
one cannot be sure that the solution is the one that corresponds to the 
infimum of the weighted sum of squares of the residuals. It is always 
possible, depenaing on the starting approximate solution and the nature 
of the particular function involved, to converge to a local minimum. 

To assure the global convergence one would have to use several and 
widely differing starting approximations and establish that the algo- 
rithm always converges to the same minimum. This can hardly eve*^ be 
accomplished in practice due to the large number of observations in- 
volved, but then again one has in most cases a very good idea of what 
values the parameters take on, and therefore in practice such dismal 
cases are scarce [Hamilton, 1964]. 

What is very real, however, is the fact that the DOC process as 
applied herein, and almost everywhere else, is neither a least squares 
adjustment nor a minimum variance one in the standard statistical sense 
of these terms [Rao, 1973]. The required partial derivatives in the 
Taylor expansion are determined in part numerically from the integration 
of the variational equation of state rather than from some well-defined 
analytical expression. It is then possible, even probable, that incor- 
rect modeling or omission of the effect of dominant forces in the orbit- 
al model will result in an incorrect or strongly biased solution. The 
fact that the solution is obtained by means of the formulas for either 
of the aforementioned statistical adjustment procedures is not enough 
justification to call the result unbiased. Since we never know all the 
mechanisms that govern the orbit of the satellite perfectly, we are not 
entitled to use the term unbiased even in the cases when all the unknown 
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parameters entering the mathematical expressions of these mechanisms 
are being determined from the observations themselves. The best that 
we can probably have is a "conditionally unbiased" estimate of these 
parameters, the condition being that the assumed model reflects reality 
to a degree that exceeds the effect of computer round-off errors in the 
final result. Obviously, with the limited knowledge of the model, it 
is to our advantage to use the available observations in ways that will 
minimize the effect of this deficiency in the resulting solution, which 
is one of the objectives of this investigation. Nevertheless, the condi- 
tional unbiasedness of tne results still holds, and it would be unrealis- 
tic to advocate the opposite just because technically the same mathemat- 
ical formulas are used in both adjustment algorithms. 

We finally want to address here an issue that some might object to, 
that is, the use of quantities (the simultaneous range-differences) as 
observations in the estimation which are not really observed, but rather 
inferred from observed quantities (the ranges). If we wanted to avoid 
this question but still be able to use the observed ranges in a differ- 
encing mode, we would have to modify our mathematical model in a way 
that the difference of simultaneous ranges and the parameters satisfy 
condition equations involving both. This mathematical mode, known as 
"combined model" [Uotila, 1967], will result in a solution that is iden- 
tical with that obtained from the simpler model where the "pseudo- 
observed" SRD's are written as a function of the parameters alone 
(observation equations model ) . 
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4.2 Differential Relations Between the Observable and the Parameters 
Let F be a vector function in X, continuously differentiable. A 
first-oroer approximation of T over a differentially small region about 
X = can be obtained by expanding F in Taylor series and retaining 
terms up to and including the first derivative: 


F(X) : F(Xj) + 


3X 


(X - X^) 


( 3 ') 


Using (91) we can linearize the SRD function in order to be ahlc to 
adjust them in the DOC process. 

Let S denote the vector of Cartesian CES referred coordinates of 
the satellite at the instant of an observation 6p, and 6i , G 2 , the cor- 
responding coordinate vectors of the observing stations in the same sys- 
tem. We can express 6p as 


[6p^l - [(S^. - G2)^(S^ - G2)]^ - [(S^ - Gi)^(S^. - Gi)]^ (92) 

where i is the number of observations. We can now identify [6pj] as 
F(X), and X as the vector containing the station and satellite positions 
along with several other parameters relating to the orbital model and 
the CIS to CES and CTS to CES transformations. 

The relationship between these frames of reference is expressed 
through the orthogonal rotations for precession (P), nutation (N), 
earth rotation ( 0 ), and polar motion (C) as 


S. = [NP]R. 

(CIS ^ CES) 

(93) 

G . = [C0]T Uj 

(CTS ^ CES) 

(94) 


Because of the similarity of the tv;o terms in (92), we need only form 
the partial derivatives for the first term; the ones for the second can 
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then be readily obtained by changing the index in the G. vector. The 
partial derivatives for 6p^ can then be formed by differencing these 
two partials. From the above the only elements of the transformations 
that will be treated as unknown are contained in the matrices C and 0. 
From [Mueller, 1969], 


C s 


and 


0 = 


1 0 X 

0 1 -y 

-X y 1 

COS0 sine 0 
-sine cose 0 

0 0 1 


(95) 


(96) 


where x,y are the coordinates of the celestial pole with respect to a 
local tangent plane coordinate system with its origin at the North CTS 
pole, its x-axis on the A=0° meridian, and its y-axis on the A=270° 
meridian, and e is the angle of rotation between the first axis of the 
CTS (origin of longitudes) and the first axis of the instantaneous CES. 
The X and y will be conside»"ed as parameters of the problem to be deter- 
mined in the adjustment. We can also include the rate of change of e 
as a parameter in hopes of determining "lerigth-of-day" variations; 
however, for reasons explained in [Van Gelder, 1978], this is best deter- 
mined from observations with alternate techniques such as VLBl (Very 
Long Baseline Interferometry), rather than SLR. 

The orbit is adjusted in terms of corrections to the initial ap- 
proximation of the satellite's orbital elements at a fixed epoch. The 
instantaneous discrepancies at the epochs of observation are related 
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to these parameters through the transition matrix Y, obtained from the 
integration of the variational equation of state as formulated in Sec- 
tion 2.2.1, equation (6). This enables us to write: 


dR. = Y[dRoi dRo]^ (97) 

— ' — T 

where [Ro I Ro] is the initial epoch satellite - state vector. 

Working with one range p.. at a time, where 

' %J 

^ [(s, - (S, - Gj)]i (98) 

and letting T. ■ denote S- - 6. for brevity, we use the chain rule of 
1 J 1 J 

differentiation to obtain: 



Using now (93) and (94) and the definitions for C, 0, and T.., the 

' J 

required partial derivatives for evaluation of (99) are determined ex- 


pl icitly : 





^ = -3[C0]I 


auj 





( 100 ) 


( 101 ) 


( 102 ) 


89 



3S 


— ' 3 [NP]; 


aR, 


ORiaWAL PAC' '■ 
OF POOR QUALITY 


L\3Rq/ I \3R()/ _ 
= -3[I]3 

36. 


= 3[Y], 


3Gj _ 


3x 

3y 


-Uj2 cose 

-Uj3 

L 


■Uj3 

Uj2 COS0 
-‘^32 


(103) 

(104) 

(105) 

(106) 


(107) 


Expressions (100) through (107) can be substituted in (99) to 
obtain the explicit first-order differential of p^. with respect to the 
parameter vector [U. I Rq ! Ro 1 x I y] . Writing the resulting equa- 

J 

tion for j=l and j=2, and subtracting the two, we obtain the differ- 
ential 


d(5p^ - dp^2 “ ^P-ji (108) 

which corresponds to the SRD of equation (92). 

Assuming now that our initial approximation for the unknown param- 
eters is close to their true values, we can use this differential ex- 
pression (108) in connection with the Taylor expansion (91), where the 
following equalities are identified: 


— _T ' -T ' -T 

X = [U^ : r! ! Ri 



(109) 


F(X) = [6p^ + e.] 


( 110 ) 
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with denoting errors to be estimated during the adjustment. Note 

that in practice there will he a number of x,y parameter pairs in a 
problem since these represent averages of the coordinates of the pole 
over a predetermined time interval (e.g., five-day averages). With 
being the initial approximation to the parameters, and lettinc 3 = X - 
X^, we get: 

6 = [d0{ i dfil i dRi i dx i dy]"^ (111) 

J 

We will use the notation A^. to indicate the ith row of the design 
matrix which consists of the partial derivatives with respect to the 
parameters X (arranged in precisely the same order). In addition, we 
let Y denote the residuals, estimates of the errors e, and d, the dis- 
crepancy vector: 

d = F(X^) - [6p.] (112) 

We can now write the linearized mathematical model in terms of the 
established notation, which results in the set of observation equations 
to be adjusted: 

[r.] = [A^-j][ej] + [d.] (113) 

The explicit expressions for the elements of the design matrix are 
further developed in Appendix C. 
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4.3 Estimable Parameters 

The concept of "estimabil ity" or "estimable parameters" was first 
introduced by R.C. Bose [1947] in an attempt to expand the applicability 
of the well-known theorem of Gauss-Markov [Rao, 1973]. In [ibid.] the 
estimabil ity of the parameters or parametric functions in an estimation 
problem is examined through the rank of the design matrix of the experi- 
ment. Because of the fact that a matrix can become "algorithmically" 
rank deficient due to ill-conditioning [Forsythe and Moler, 1967], there 
has been some confusion in the past as to the status of geodetic param- 
eters estimated from dynamic space techniques. 

The truth of the matter is that the status of a parameter in this 
respect is determined by the underlying physical reality and not by the 
numerical entries of a matrix or their interrelationship. The scienti- 
fic interest in this subject is reflected by some rather extensive 
literature, the most recent ones being [Van Gelder, 1978, Grafarend and 
Livieratos, 1978; Grafarend and Heinz, 1978; Pavlis, 1979]. None of 
these investigations, though, has looked at this issue from the physical 
point of view. The estimability of a parameter should be determined by 
two simple factors. Eitner the observations contain information about 
the parameter, or the model contains information, or both. 

It is a trivial exercise tc show that »^ange measurements and linear 
measurements in general contain only scale information and are indeoen- 
dent of the coordinate system that is used in the parametnzation of 
the problem. On the other hand, the physical model describing the dy- 
namics of the satellite contains information that is enough to define a 
particular coordinate system. If the harmonic coefficients for the 
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geopotential are assumed known, then the satellite orbit becomes 
sensitive to the coordinate system in which these coefficients are 
referenced. It is only a peculiarity of the physical figure of the 
Earth that two out of the th»ee principal moments of inertia are nearly 
equal [Heiskanen and Moritz, 1967], namely, the equatorial moments; 
the definition of the origin of longitudes, therefore, is a very weak 
one. This is, in fact, the reason why in dynamic solutions with only 
metric measurements the longitude of one of the participating stations 
is kept fixed. This is not the only solution to overcome a case of pure 
ill-conditioning, but it is the most popular and the simplest to apply. 

It is thus obvious that there is no rank deficiency in the dynamic 
problem of satellite geodesy conditional on the fact of a finite expres- 
sion for the geopotential, but rather an extreme ill-condition due to 
the aforementioned reasons. With that in mind we can further investigate 
the interaction between parameters of interest to determine which of 
them are separable in a simultaneous adjustment. This can be best 
accomplished by examining the information required for their determina- 
tion. 

4.3.1 Information required for the determination of the problem 
parameters. 

The goal of this investigation is primarily the determination of 
interstation baseline distances and variations in the coordinates of 
the celestial pole. The formier are obtained from the estimated station 
coordinates, while the latter are determined as additional rigid body 
rotations of the station network with repsect to the satellite orbit, 
in addition to the modeled rotations included in the transformation 
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between the CTS and the CIS systems. The estimation will be based on 
simultaneous range-differences reduced in a long-arc dynamic mode solu- 
tion. We examine here the complete observation equation with parameters, 
the station positions, the coordinates of the pole, and the initial state 
vector of the satellite to determine which of them are separable and 
therefore design experiments in which those parameters of primary 
interest will be estimable. 

For the sake of brevity, the form of the observation equations 
used in this section is the "body-fixed" form presented in Appendix C. 

In the following, the prefix l will denote differences of coordinates, ? 

W will denote the variational partials matrix rotated in the CTS system 
with W(j) denoting the jth column of that matrix. Subscripts will refer 
to station positions and superscripts to satellite positions. The 
letters X, Y and Z will denote CTS coordinates, and coordinate differ- 
ences are always taken between a satellite position and that of the 
observing station, i.e., 

aXJ = X^ • X 2 (113) 

We write here two observation equations, each from a different 
pair of stations to a different satellite position. The parameter list 
is in the following order: 

(1) Station coordinates X,Y,Z for each of the stations, 

(2) Satellite state vector at initial epoch, and 

(3) X, y coordinates of the pole. 

With three parameters for each of the four stations, six for the state 
vector, and two for the pole, there are twenty parameters in total, and 
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the design matrix A given in equation (114) will have twenty columns. 

In identifying linear dependencies between the parameters, we will use 
the notation A(j) for the column of A. We can separate A into three 
submatrices according to the three major groups of parameters as pre- 
viously stated: 

A = [Ag i A3 i Ap] (115) 

with Ag consisting of the partial s with respect to ground station coor- 
dinates, A3, partials with respect to the satellite state vector, and 
Ap, partials with respect to the coordinates of the pole. From the 
definition of these groups and (114), we can state that there are no 
linear dependencies betwe' columns of the design matrix within the 
same group. We concentrate, therefore, in finding the dependencies, if 
any, between columns of the design matrices ^mong different groups. We 
are investigating here whether it is possible to find a set of constants 
{Cl, .... C2o)> not all zero, such that the following equality holds: 

CiA(l) + CzA(2) + ... + C2 oA(20) = 0 (116) 

At first glance it would seem as if the columns of the A3 sub- 
matrix can be written as combinations of those of Ag. This is not so, 
though, since the W matrix is different from one observation to the 
next; and therefore its elements w^.j which are used in for :ing A3 differ 
too (which is obvious by the different superscripts in (114)). In the 
case of Ap and Ag, however, one can easily write the following relation- 
ships between their columns for the first observation: 

(-Z;)Ag(4) 1- (X2)A|,(6) + (-Zi)Ag(l) + (Xi)Ag(3) Ap(l,l) (117) 

(-YJA|,(6) + (Zj)Ag{5) + (-Yj;Aj( 3) + (ZJAg(2) = Ap(l,2) (118) 
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and for the second: 

(-Z4)Ag(10) + (X4)Ag(12) + (-Z3)Ag(7) + (X3)Ag(9) = Ap(2.1) (119) 

(-Y4)Ag(12) + (Z4)Ag(ll) + (-Y3)Ag(9) + (Z3)Ag{8) = Ap(2.2) (120) 

Observing now the location of zero elements in Ag, we can combine the 
above in one equation that shows the existing sought for set of constants 
C to be: 

C = [-Zi Zi -Yi+Xi -Z 2 Z 2 -Y 2 +X 2 -Z 3 Z 3 -Y 3 +X 3 -Z 4 Z 4 -Y 4 +X 4 

0 0 0 0 0 0 -1 - 1 ] ( 121 ) 
so that: 

aJ = [4,] (122) 

The zero elements of C correspond to the columns of A which constitute 
the submatrix A^. Dependencies between A^ and Ap do not ex ‘it, for if 
they did we would come to the contradiction that Ap which is a linear 
transformation of Ag can be written as a transformation of Ag and at the 
same time A^ can be independent of Ag. 

We come to the conclusion, therefore, that we cannot separate the 
station parameters and the ones for the pole in a simultaneous adjust- 
ment. Theoretically, this dependence would be broken if (122) did not 
hold even for just one row of A. If that is the case, we assumed that 
x,y are known for the first (say) observation, and we set the correspond- 
ing partials A(l, 19), A(l, 20) equal to zero. The problem is that if 
we were to check the column independence of the resulting A matrix 
using (122), we would find that the result differs from zero only 
slightly. That means the corresponding parameters have extremely strong 
correlations, and the normal equations will be algorithmically singular. 
If the above is repeated for several rows of A, the condition of the 
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normal equations is improved, but one must realize that in this case the 
station coordinates are being determined from these first observations, 
the rest of them contributing very little due to the confounding between 
the two groups of parameters. 

Conceivably, we could assume that the coordinates of the pole are 
known from previous solutioi.s, for a sufficiently long period of time 
so that data available over that period can be used for the estimation 
of station positions alone. With our present capabilities in laser 
ranging, it would be required that data over more than a month's inter- 
val be used for this purpose, in order to achieve sufficiently accurate 
station positions for subsequent estimation of the global motions of this 
network. Objections that one can raise against this practice is the 
bulk of computations that need to be repeated in every solution and the 
fact that the coordinates of the pole estimated in each of these solu- 
tions refers to the CIS defined by the simultaneously estimated station 
positions. That is to say, every solution defines a new CTS. This is 
obviously the most undesirable of the two side effects of a simultaneous 
solution. It would therefore be more meaningful to do a separate solu- 
tion for the stations from that for the coordinates of the pole. Adopt- 
ing the resulting CTS at some epoch, we could then monitor the rotations 
of the defining network of stations with respect to the celestial pole 
by 'eans of independent solutions in which the station positions are not 
allowed to adjust. 
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4.^1 Minimization of Model Biases by Use of SRD Observations 

We have already seen that the adjustment of satellite ranging data 
is based on a model that involves hundreds or thousands of parameters, 
none of which is perfectly known. Moreover, it is practically impossi- 
ble to include all of them as unknown parameters of the problem. We 
therefore resort to the option of adopting a fixed value for a number of 
them. Such a practice obviously biases the results of the adjustment 
since the adopted values differ from the ever unknown true values of the 
fixed parameters. From the description of the ranging model it should 
be clear that the most parameters which are held fixed are those involved 
in the determination of the satellite orbit (e.g., geopotential coeffi- 
cients). Their errors affect the quality of the orbit directly, and 
when the model value for the range is computed from this orbit to com- 
pare it to the observed one, the errors propagate into the discrepancy 
vector and thereby in the solved for parameters. We will see now how 
these biases can be diminished if we take advantage of the simultaneity 
of the observations and use them in the simultaneous range-difference 
mode. 


Assume that we have two sets of range data in which each observa- 
tion from the first is taken at the same time to the same satellite posi- 
tion as the corresponding one from the second set. If we were to differ- 
ence these data sets we would obtain the SRD data as we have described 
them in Chapter 3. The linearized observation equations for these two 
data sets can be written as 
ri = Ai 6 + di 


(123) 


^2 ~ A 23 d2 
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and the solution for the parameters 6 is 

3 “ “(a 1 Ai +aJ A 2 ) ^(Ai J^i^di + A 2 ^ 2 ) (124) 

We have considered here the case of common parameters 3 for both data 
sets since the data are taken on the same satellite arc, and as we 
mentioned above the major source of bias is the computed satellite 
orbit. 

The discrepancy vrctor d can be written out in terms of its 
components as 

d, = F.(b^) - (125) 

where (^^-(3^) is the model computed range and p° the one observed. If the 
fixed parameters were fixed at their true values, then the computed range 
would have a different value, p^. That would, of course, be the desired 
value also, although this is practically impossible. The term ^^^(3 q) 
then has two components, one being p^ and the other being the bias b^. : 

F,(6o) = Pj + b. (126) 

Using (126) in (124) and denoting with N the matrix of normal 
equations, the resulting expression for the solved for parameters 3 
becomes 

3 = -{N'^aI I?(p9 - p°)i + aI - P°)2] + 

+ N'^[aI zI^ bi + aI z? b2] (127) 

The first term in (127) represents the proper adjustment in the 
parameters 3 based on the information in the observations, while the 
second is the bias term in 3 due to errors in the adopted values of the 
nonadjusting parameters of the model. A secondary and much less serious 
effect of the erroneous model parameters is the error in the elements 
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of the design matrices and Aa which in general are functions of these 
parameters. These errors, however, are more important in computing the 
covariance matrix of the adjusted parameters than in estimating e. 

We will now derive the corresponding adjustment equations for the 
SRD mode. Since the order in which the observations have been arranged 
is not important in the adjustment, we may assume that the two data sets 
are already in correspondence, i.e., 

Pli P 2 ^- in the sense that: 6p^. - P 2 ^ - P^^ (128) 

From (123) and (128) then it follows that the linearized observation 
equations for the SRD mode are: 

Aa + d = r (129) 

with 

A = [A2 - Ai], 

d = d 2 - di, and (130) 

r = r 2 - ri 

The minimum variance adjustment estimate for the solved for parameters 
will then be 

3 = -{[(aI- aI)(Zi+ E2)’M^2 - Ai)]'^ (aI-aI)(Ii+ l2)"Md2-di)} (131) 

where we have used Zi + z,, ^s the variance-covariance matrix of the SRD 
pseudo-observations as obtained through error propagation on the basis 
of (128). 

Writing the discrepancy vector in its components again, we obtain 
the following: 

d s d2-di = (p5-p?) - (p?-p?) + (b2-bi) = 6p^ - 6p° + (b2-bi) (132) 

Denoting by N the matrix of normal equations in (131) and considering 
(132) we can write 
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0 = -{N'H(aI- aI)(Zi+ E2)‘M6p^- <5p°)] 

+ N"^[(aI- aI)(Ii + Z2)'Ubz - bi)]} (133) 

for the estimate of the parameters' adjustment from SRD observations. 

The last term in (133) is again the bias term, only this time it involves 
the difference of the bias in each of the computed ranges. Naturally, 
the two bias components bi and bz do not always have the same sign or 
magnitude, so we cannot outrightly set the last term to zero. We will 
examine though the behavior of this difference in comparison to its indi- 
vidual components by means of a simulation study. 

4.4.1 Simulation study for bias propagation characteristics. 

The computation of 6p^ is based on the coordinates of the observ- 
ing stations and those of the satellite at the epoch of the observation. 
We will introduce known biases in the satellite coordinates and then 
examine how these biases and to what extent they affect the computed 
ranges Pi and pz, as well as the corresponding SRD's 6p^. To simplify 
the computations, we will assume a spherical earth model and a satellite 
orbiting at the mean altitude of Lageos on a circular orbit. These 
assumptions are well justified since they are not too far from the real 
situation, and we are only after order of magnitude of the biases rather 
than exact numbers. We have computed the biases on the intersection 
points of a l''xi° grid covering the area around the observing stations 
and then plotted the results to ease their evaluation. It is common to 
state orbital errors in three directions, the radial, the along-track 
and the across-track, but we have chosen to use the radial, latitudinal, 
and longitudinal directions. This simplifies the computations without 
altering the results, and it also has the advantage that the bias-surface 
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plots are applicable for a much wider class of orbits, while otherwise 
we would have to restrict ourselves to those having the same inclination 
with Lageos as well as the same mean orbital altitude. Although we have 
tested several station configurations in terms of absolute position on 
the earth's surface, as well as relative to each other, we will present 
and discuss here only two cases which are representative of the complete 
set. In case A the observing stations are 2000 km apart (chord distance 
between them), while in case B we decreased the distance to 200 km. For 
both cases we have taken the first station to be at latitude 40° N and 
longitude 0°, while the second occupies four diffcirent positions so that 
the azimuth of the great circle arc connecting the two is 0°, 30°, 60°, 
and 90° counting clockwise positive from the meridian of the first station. 

The satellite coordinates at the grid points are computed from the 
spherical coordinates r, cji, X using 


X = r coS(j) cosA 

Y = r COS4) sinX (134) 

Z = r sin;;/ 

and the biases in X,Y,Z are obtained from the adopted values in the r, 

4>, X system using the following transformation which follows from differ- 


entiation of (134) : 



where we have used the following substitutions 


(135) 
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sin4) = 

Z/r 

cos(f) = p/r 

(136) 

sinX = 

V/P 

cosX = X/p 



The values of Ar, Ac|), and AX used in the simulations are 
Ar = 2.00 m 

A4> = O'.'Ol 
AX = -0':02 

These values were arbitrarily chosen; their order of magnitude 
though reflects the present level of stability in the satellite orbit 
defined and maintained reference frames. 

4.4.2 Analysis of the simulation results. 

Using equation (135) we have separated the biases into their com- 
ponents AX , AX,, AX„, etc. so that we can study individually the propa- 

(f) A I 

gation characteristics c*' each one of them and their effects on the com- 
puted ranges and SRD's. 

This arrangement resulted in three sets of biases (the radial, the 
latitudinal, and the longitudinal) for each station configuration consid- 
ered. These biases are, of course, three-dimensional functions that 
depend on the absolute as well as the relative positions of the observing 
stations and the observed satellite points. An optimal way of displaying 
their features and characteristics is to create contour plots in the 
regions of interest. This is what is presented and discussed in this 
section. The contour plots which are discussed here are based on the 
results of the simulations described in Section 4.4.1. 

The quantity which is plotted is the bias in the range/SRD due to 
the biases introduced in the orbit. In all cases this bias is plotted 
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in centimeters. The following example will clarify the use of these 
"bias surface" contour plots. 

From Fig. 7(a) we find that when station 1 is ranging to Lageos 
at a point with subsatellite coordinates ({>= 37° N and X= 10° E, the 
bias in the computed range due to the 2.00 m radial bias in the orbit 
used will be 1.98 m or 198 cm. Fig. 8(a) shows that when the same satel- 
lite point is observed from station 2, the bias in this case is about 
1.88 m. When the two ranges though are subtracted to create an SRD ob- 
servation, the resulting bias is only 10 cm (!), as can be verified from 
Fig. 11(a). 

Inspection of the bias surface contour plots leads to a number of 
interesting remarks. The radial bias surface is bell shaped with its 
apex on the observing station's radius. The latitude and longitude bias 
surfaces exhibit an antisymmetry, the former with respect to a line of 
almost constant latitude (equal to that of the observing station), and 
the latter with respect to the station's meridian. The form o^ these 
three surfaces depends on the absolute position of the station only in 
the case of the latter two, and in this case the one for latitude depends 
only on the latitude of the station. In any event, their shape changes 
very slowly and since a coobserving pair of stations should not be more 
than about 2000 km apart (in case of Lageos), for all practical purposes 
we can assume the corresponding surfaces to be the same. Finally, as 
far as the smoothness of these surfaces is of concern here, we can see 
from their contour plots that at least in the vicinity of the observing 
stations the radial bias surface is by far the most flat of the three, 
the other two exhibiting considerably stronger gradients. 
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The bias surfaces for the range differences are computed exactly 
as the differences of the true SRD values from the biased ones. They 
represent therefore the difference surface of the corresponding range 
bias surfaces. On the basis of this observation and the previous discus- 
sion on the propagation characteristics of the various biases, it should 
come as no surprise that radial biases are almost completely eliminated 
when we use the SRD mode. Compare, for instance, the contour levels 
between Figs. 7(a) and 11(a). In addition to this we also note that the 
level to which this bias can be eliminated depends jn the distance be- 
tween the two stations as well as the relative location of the satellite 
track and the interstation baseline. The closer the stations, the smaller 
the remaining bias in the SRD's, as one can verify from Figs. 11 and 12. 
Since the two surfaces are nearly the same in the vicinity of the sta- 
tions, their difference will be smallest in the area between the two 
stations and close to their baseline. If one now considers the fact that 
as the interstation distance increases, the area in which simultaneous 
observations are possible "shrinks" towards the point amid the two sta- 
tions, it becomes obvious that the SRD mode has a clear advantage over 
simple ranging in the case of radial bias in the satellite orbit. 

The situation is quite the opposite in the base of latitude and 
longitude biases. Since the range bias surfaces in this case have dif- 
ferent signs in different areas, it is possible that for some areas the 
biaseswill increase in the SRD mode while they will still decrease in 
others. This is indeed what happens in the area in between the two sta- 
tions as Figs. 7(b), 9 and 13 show for the latitudinal bias, and Figs. 
7(c), 10 and 15, for the longitudinal case. If, however, the length of 
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the interstation distance is decreased so that most of the coobservable 
satellite positions lie outside the critical area, then the SRD mode 
will again be biased to a much more limited extent than the ranging mode. 
This is illustrated by the 20C km baseline example in Figs. 14 and 16 
which should be compared to the range bias surfaces shown in Figs. 7(b) 
and 7(c). We do not have to go to the extreme of using only very short 
baselines since the high altitude of Lageos allows for a rather extensive 
area of coobservable satellite positions even when the interstation dis- 
tance is quite long. With a spherical approximation, the radius of 
Lageos observability around a station is about 60°. That means an area 
that is almost one-third of the total of the globe. Hence even with 
longer baselines, the SRD mode can be considerably less affected by 
biases in the latitude and longitude directions compared to the ranging 
mode if the data are collected in a region that excludes the immediate 
vicinity of the baseline. 

In general, all three types of biases will affect the computed 
range differences, and we thus have to select our data in a way that all 
of them will be minimized simultaneously. Since the radial bias is mini- 
mized independent of the satellite position, the guidelines are set by 
the requirements we mentioned above for the minimization of the latitud- 
inal and longitudinal biases. 
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Fig. 7 Range bias surfaces for the fixed station 1. 
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Fig. 13 Latitudinal bias surfaces for SRD's from four 2000 km baselines 
(values in centimeters). 
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Latitudinal bias surfaces for SRD's from four 200 km baselines 
(values in centimeters). 
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Fig. 15 Longitudinal bias surfaces for SRD's from four 2000 km baselines 
(values in centimeters). 
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4.5 Sensitivity Analysis of SRD Observations 

4.5.1 Introduction. 

In designing an experiment the minimization of bias effects with- 
out overparametrizing the problem is an important factor, but it is 
neither "the" major factor nor the only one to be considered. A more 
important and challenging issue is the selection of the proper data that 
would yield the most accurate parameter estimates with the least amount 
of computations. Not all of the available observations will contain the 
same amount of information about the parameters in general. That does 
not mean that we should not use these observations, but if the computa- 
tions are tedious and expensive in terms of required computer time, then 
we might want to weigh the benefits from the use of these observations 
against the increased computation effort and cost. 

Optimization of a design in classical geodetic networks has been a 
popular topic with a very rich literature. In the case of space syst&n 
networks, though, the problem becomes extremely complex and practically 
intractable on an observation-by-observation basis due to the possible 
thousands or even hundreds of thousands of observations connecting the 
space subset and ground subset of network vertices. We have to resort, 
therefore, to a geometrical analysis of "categories" of observations 
rather than a one-by-one examination. Since the observations are col- 
lected from individual satellite passes over the network stations, it 
seems natural to use the segregated observations over each pass as an 
entity whose optimal position relative to the network is sought. Opti- 
mality again is judged by the contribution of information about the prob 
lem parameters. Going from individual observations to individual passes 
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reduces the amount of work tremendously. Further savings can be 
achieved by classifying the passes according to their direction (e.g., 
North-South, East-West, etc.) and the maximum elevation that the satel- 
lite reaches with respect to the station's zenith (e.g., overhead passes 
horizon passes, etc.). We therefore end up having to deal with only a 
small number of categories of obse’^vations which is a much more manage- 
able problem than what we originally started with. 

4.5.2 Optimal designs for baseline estimation. 

Our primary interest here is to find optimal estimates of the in- 
terstation chord distances (baselines). An extensive investigation for 
the range observable has been published in [Pavlis, 1979]. Since ranges 
and SRD's contain the same type of information, the conclusions of that 
study are valid for our problem too. Some of the most important results 
are the constancy of the relative precision of the system over a wide 
range of baseline lengths, the accuracy dependence of the estimated base 
line on the orientation of the pass relative to its direction, and the 
independence of that estimate from biases in the adjustment's parameter 
estimates. 

The first result is important from the point of view that we are 
not restricted to use baselines of equal length throughout our design. 
The last one is of importance also, especially in cases where due to in- 
sufficient data or unfavorable distribution a desired unbiased estimate 
for all or some of the parameters cannot be found. The usual practice 
in this case is to apply prior information (if available) and use the 
Bayesian adjustment process, whereby we deflate the variance of some of 
the parameters in exchange for the unbiasedness of the solution [Pavlis, 
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1979]. In our study, for instance, we have already see». that the statio'" 
position and the coordinates of the pole are inseparable parameters. 

When the one is being estimated, the other mu';! be kept constrained to 
already known val js, and i^ is, of course, desired that the estimates 
be least affected by any bias introduced by the constraints. 

The relative orientation between baselines and satellite passes 
will be reexamined here, since the peculiarities of the SRD observations 
shed more light on the relationship between observables and parameters. 
The results of the aforementioned study are still valid, that is, satel- 
lite passes parallel to the estimated baseline should be preferred to 
those which cross it at almost a right angle. Since the SRD's, though, 
are the differences of ranges emanating from the endpoints of the base- 
line and directed to the same satellite position, they can be expressed 
as functions of the baseline length directly. From Fig. 17 and some 
elementary geometry we see that if d is the length of baseline 1-2, and 
6 is the observed range difference S1-S2, then 

AB = d cosa) (137) 

where from the triangle ISD we find 

(i) = a + ij; (138) 

From the triangle ICE we obtain 

AB = CE = Cl cosij; (139) 

and since 

SC = S2 => Cl = SI - S2 or Cl = 6 (140) 

Therefore 


AB = 6 cosijj (141) 

which, upon substitution in (137), yields 
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f^ig. 17 Satellite pass - baseline geometry for simultaneous 
range- differences. 
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(142) 


. COS(g 
COSi|; 

It becomes obvious from (142) now that if the observations were collected 
near the satellite point equidi:>tant from stations 1 and 2, then their 
values woula be near .zero and the determination of d very poor. On the 
other hand, as a -*■ 0° or 180', which happens as S -*> Ni or Nj respectively, 
then 6 d and uncertainties in the values of a and tp have no effect on 
the determination of d. In thiscasewe have assumed that the baseline 
lies on the plane of the pass; therefore its direction intersects the 
satellite trajectory at Ni and Nj. Although this special case illustrates 
to satisfaction the relationship between observables and parameters of 
interest, it is hardly ever possible for it to happen in practice; and 
even i ' it did, we would have no means of knowing a priori anyway. Ttie 
more realistic case is that of pass number two in Fig. 17, where the 
orbital plane and the direction of the baseline are nearly parallel. In 
this case the points Nj and Nj do not exist, but the observations which 
are collected at low elevation angles, near the points where the horizons 
of the two stations intersect with the satellite trajectory, will still 
be the ones with the greatest amount of information about the baseline. 

For stations which are not too far apart (baselines of a few hundred 
kilometers), the laser rays will travel through almost the same atmospher- 
ic layers and any errors from an incomplete refraction model which is 
always a limiting factor in low elevations, will be highly correlated 
and therefore cancel out in the differencing. Using the SRD mode for the 
reduction of laser ranging data, therefore, we can make use of the low 
elevation observations that are normally edited from the range adjustment 


122 



and at the same time enjoy the minimization of biases due to errors in 
the mathematical model which this approach offers. 

4.5.3 Op:imal design for the estimation of the motion of the pole. 

Basel in-<^ die, of course, not the only parameters of interest. 

The coordinates of the celestial pole are equally important here as well. 
The former are parametric functions of the adjustment parameters, and as 
such we could not use the sensitivity matrix of the design to deduce the 
optimal experiment for their estimation. For the latter, though, this 
approach can be taken since they enter the adjustment directly. The 
goal of this investigation is to find an experimental design for SRD mea- 
surements which will result in a set of normal equations with an associ- 
ated matrix being as close to a diagonal matrix as possible and with as 
large diagonal elements as possible. Such a normal matrix will, of course, 
result in small parameter variances and insignificant correlations among 
them. To put it in fewer words: an orthogonal design. In the case of 

a truly orthogonal design, the columns of the design matrix A are orthog- 
onal, and since the nomal matrix is the Gramian of those vectors, the off- 
diagonal elements are all zero. When exact orthogonality cannot be met, 
we must try to stay as close to such a design as possible. 

One way to do this, and probably the most illustrative, is to exam- 
ine the variations in the sensitivity of the observable with respect to 
the parameters as a function of its various possible realizations. In 
other words, compare the elements of the des’qn matrix for each of the 
parameters at all possible observation points. Since each SRD observa- 
tion involves two station locations and one satellite position, there are 
an infinite number of variations of this three-point configuration to be 
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examined. We have again a practically intractable problem. The pole 
coordinates, though, are basically two orthogonal rotations about the x- 
and y-axes, and one can therefore use an intuitive approach to identify 
a small number of configuration variations that would be enough to sup- 
port the analysis. Based on the form of the sensitivity matrix elements 
that correspond to the (x,y) parameters and the definition of the coordi- 
nate system in which these parameters are referenced, we conclude that 
two absolute locations for the observing stations that could give infor- 
mation about our problem are those near the meridional planes on which 
the X and y axes lie, i.e., the X = 0° and X = -90° meridians. 

The rotations for the motion of the pole are about axes that lie on 
the equatorial plane, and therefore their effect on the coordinates of 
stations increases as the latitude of these stations increases. From 
the practical point of view though, we cannot expect to have laser sta- 
tions in near polar latitudes, and we have therefore limited our tests 
to areas where most of the currently operational stations are located. 
This study follows very closely the setup for the bias propagation study. 
We have again used spherical approximations and a circular orbit at 
Lageos' mean altitude, and we have plotted the values of the sensitivity 
coefficients for each of the parameters at the intersection points of a 
l°xl° grid surrounding the coobse»^ving stations. In order to examine the 
dependence of the observable's sensitivity on the relative positioning of 
the stations in the case of the SRD mode, we have used various baseline 
lengths and azimuths to determine the second station with respect to a 
fixed position of the first one. From the analytical expressions for the 
sensitivity coefficients, we can gather that their numerical values will 
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in general increase as the coordinate separation between the stations 
increases. This observation has been verified by the numerical tests 
for different baseline lengths. It has also been observed that the abso- 
lute maximum of these sensitivity curves is encountered near the midpoint 
of the baseline, and their gradient is minimum in a direction nearly per- 
pendicular to the baseline at its midpoint. These observations and the 
fact that the simultaneity requirements constrain the actual baseline 
lengths to not more than about 2000 km indicate that a 1000 km separation 
would be the optimum for obtaining a sufficient number of observations 
with high enough sensitivity for a precise parameter determination. We 
have therefore included here only the sensitivity plots that refer to 
this particular case. 

The sensitivity surfaces are shown for the range observable from a 
station at 40° N latitude in Fig. 18 for 0° longitude, and Fig. 19 for 
X = -90°. It can be readily verified from these plots that as the sta- 
tion is moved in longitude by 90° to the west of its original position, 
the sensitivity surfaces undergo a 90° counter-clockwise rotation about 
the station's geocentric radius. The end result is that the sensitivity 
surfaces for x and y at 0° longitude are identical to those of y and x 
at 90° W longitude except for a sign change in the coefficients for x at 
the second location and those for y in the first. 

Inspection of other cases where the stations are located in between 
the above two meridians shows that indeed the shift in sensitivity from 
one parameter to the other is a smooth operation and on the A = 45° merid 
ian (or A = -45°) the range observable is equally sensitive to both param 
eters. It is also worth noting here that the sensitivity of tiie system 
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Fig. 18 Sensitivity surfaces for the coordinates of the pole 
in the case of a station at X = 0°, c}) = 40°. 
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(a) x-coordinate 



(b) y-coordinate 


Fig. 19 


Sensitivity surfaces for the coordinates of the pole 
in the case of a station at X = -90°, 41 = 40°. 
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is in general increased with respect to either of the parameters as the 
observed satellite positions are selected farther away from the station's 
zenith. In the case of ranging, therefore, we are in the unfortunate 
situation that we have to use very low elevation observations if we want 
to increase the sensitivity of our system with respect to the motion of 
the pole. As it has already been mentioned, these observations have the 
largest uncertainties due to incomplete modeling of the atmospheric re- 
fraction effects and since the ranging mode cannot eliminate biases as 
the SRD mode does, such poor quality observations are hardly ever included 
in the range adjustment. 

These last few observations, when first noted in the initial steps 
of this investigation, prompted us to examine the potential use of the 
SRD mode for the estimation of the coordinates of the pole. Based on the 
sensitivity surfaces of the previously used station pairs we computed by 
differencing the sensitivity surfaces for the SRD mode. From Fio^, 20 
and 22 we can see that for a stronger determination of the x coordinate 
the station pair must be in the vicinity of the X = 0° meridian (or X = 
180°), and of all possible baseline orientations, that nearest the North- 
South direction will result in the highest sensitivity possible. Simi- 
larly, Figs. 21 and 23 indicate that the strongest determination of the 
y coordinate will result from the observations at a station pair near the 
X = -90° meridian (or X = 90°), and again the highest possible sensitivi- 
ty will be achieved if the orientation of the baseline is in the general 
North-South direction. It can be verified also that baselines in the 
East-West direction will contribute to the sensitivity of the model too, 
but in a reverse manner from that of the North-South baselines. That is 
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(a) 0° azimuth 



(b) 45° azimuth 



(c) 90° azimuth 


20 Sensitivity surfaces for the x-coordinate of the 
pole for three 1000 km baselines near X s 0°. 
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(c) 90° azimuth 


Fig. 21 Sensitivity surfaces for the y-coordinate of the 
pole for three 1000 km baselines near A s 0°. 
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(a) O'" azimuth 




(c) 90° azimuth 


Fig. 22 Sensitivity surfaces for the x-coordinate of the 
pole for three 1000 km baselines near X = -90°. 
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(a) 0° azimuth 



(b) 45° azimuth 



(c) 90° azimuth 


Fig. 23 Sensitivity surfaces for the y-coordinate of the 
pole for three 1000 km baselines near X s -90°. 
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to say, an E-W baseline near the A = 0® or 180® meridian will be sensi- 
tive to the y coordinate although almost completely insensitive to the 
X, as Figs. 21(c) and 20(c) show respectively. In parallel, for an E-W 
baseline near X = 90® or -90°, the observable is sensitive to the x 
coordinate but not to y, as Figs. 22(c) and 23(c) indicate. 

The conclusions which can be drawn from this study are rather obvi- 
ous by now. The optimal SRD network for monitoring the motion of the 
pole should consist of perpendicular baseline pairs located near the two 
meridians on which the x and y coordinate axes lie. Some savings in the 
number of dedicated laser equipment can be achieved by giving these sub- 
networks an L-shape, one station being common for both the N-S as well as 
the E-W baselines. Considering now that laser ranging is a weather- 
dependent system and the fact that with a single Lageos target there is 
a six hour gap in every 24-hour period during which a typical station at 
about 40® latitude will not be able to observe due to the earth's iner- 
tial rotation, it is only logical to plan for far more stations than the 
mere minimum. If the additional requirement of a uniform distribution of 
the observations over the globe is considered, then we should also in- 
clude subnetworks located in the Southern Hemisphere. In this manner we 
not only increase the chances of observing the target within a given time 
interval, but we can also minimize orbital biases coming mainly from an 
inaccurate geopotential model. The importance of having a uniform dis- 
tribution of data in time will be further discussed in connection with 
the operational estimation technique for the coordinates of the pole. 

What should be considered though as the major advantage of this network 
configuration over simple ranging is the fact that the resulting design 
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Is nearly orthogonal over the entire area surrounding each baseline. In 
other words, observations from each of the baselines In an L-paIr are 
only sensitve to one of the parameters, either x or y, and estimates of 
these parameters will have nearly zero correlation as the normal equations 
matrix Is almost diagonal. This cannot occur In the case of ranging 
except In very restricted designs where the Sdcelllte pass fellows either 
of the two directions on which the sensitivity with respect to one of the 
parameters Is nearly zero while 1t Is maximum for the other one (see Fig. 
18). 


4.6 Operational Approach for Parameter Estimation 

4.6.1 Estimation of baseline lengths. 

The estimation of baseline lengths In a dynamic solution Is based 
on the estimated coordinates of the observing stations. We have shown 
already that In a long-arc problem where the satellite Is tracked over 
several revolutions, the coordinates of the tracking stations are separa- 
ble from the Initial conditions for the orbital Integration (initial sat- 
ellite state vector) except for an ill -conditioning in longitude caused 
by the peculnrities of the earth's mass distribution. 

The first option that one has then is to adjust all the data col- 
lected from all stations observing the satellite in one batch adjustment 
with one station's longitude constrained to an adopted value. Such a 
solution, of course, will have to utilize a very detailed orbital model, 
since any unmodeled sources will progressively affect the solution more 
and more as the length of the arc becomes longer and longer. As the 
model gets more complicated, the computations become more complex and 



time consuming; therefore, the cost per estimated baseline Increases. 

In addition to the above, a requirement of the long-arc solution is the 
uniform distribution of the observations over the entire arc, and for 
this to be achieved the stations must have a global distribution. That 
presupposes the existence (and continuous maintenance) of a global net- 
work of stations. For various reasons, the land-sea distribution being 
the most obvious, the existence of such a network cannot always be 
guaranteed. 

An alternate approach which circumvents the aforementioned problems 
is the semi -dynamic solution where each satellite pass over a station of 
interest is treated as an independent arc. In this case the length of 
the arc hardly ever exceeds one-third of a complete revolution, and the 
orbital model therefore cannot furnish the required information about the 
origin and orientation of the underlying coordinate system. The shorter 
the arc, the less the contributed information by the orbit. Such solu- 
tions have been common practice with Doppler system users, especially 
individual ones who have no way of observing the satellites on a global 
scale. The usual remedy to the problem is the constraint of the satel- 
lite orbit over each pass to some fairly accurate known values. By doing 
so, one hopes that the introduced bias in the recovered positions of 
nearby stations will be highly correlated and cancel out during the inter- 
station distance computation. This, however, is neither guaranteed nor 
does it mean that the resulting baseline estimate is bias-free by defi- 
nition. The advantages of this method though are quite obvious. There 
are no requirements for an extensive network and the orbital arcs being 
very short, the orbital model can be simplified tremendously, since long 
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period and secular perturbations In the model have no time to build up 
and corrupt the solution. As a result, the cost of this mode of solution 
compared to the previous one (for the same amount of data) Is reduced 
substantially, In some cases by an o. der of magnitude. 

4.6.2 Estimation of the coordinates of the pole. 

The fact that the coordinates of the pole cannot be separated frum 
those of the observing stations In a simultaneous adjustment has been 
shown already in Section 4.3.1. At that point some options were examined 
for the separation of these parameters. Because of the arguments we 
raised at that point, the second option seems to be the only one which 
can produce consistent results over a long period of time. What is even 
r *e Important Is the fact that this method would result In a set of pole 
coordinates that Is compatible In sense with those currently obtained 
and published by the International organizations such as BIH. 

Since the coordinates of the observing stations are fixed to some 
adopted values (which define the CTS), errors in these values will affect 
the resulting polar motion record. As long as we always use the same set 
of coordinates the effect of their errors on the polar motion record is 
more or less the same, and we need not concern ourselves with this item 
any further. The effort should be diverted in establishing the connec- 
tion between the previously adopted CTS and the one used in the new tech- 
niques. In fact this has been one of the most pressing issues in recent 
years, and there are several suggestions on how to deal with it, the most 
recent and, in our opinion, the most straightforward appears in [Mueller 
et al., 1982]. 
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The idea of monitoring the motion of the pole by two rigid body 
rotations of the station polyhedron with respect to the CIS defined by 
the satellite orbit seems to be in agreement with the current trends and 
requirements of the interested parties. There are two items hich are 
of concern in connection with this mode of estimation--the dynamic na- 
ture of the parameters and the stability of the satellite-defined CIS. 
Since we do not have a rigorous mathematical model to describe polar mo- 
tion, we must resort to a discretization of the problem for the estima- 
tion of the state of the process. The quality of the final result will 
depend on the distribution of uhese individual estimates in time. On the 
other hand, the distribution of these estimates in time is associated 
directly with the spectrum of the process and the capabilities of the 
measuring system. The already long record of observations has establish- 
ed the fact that there are two dominant frequencies in the spectrum, an 
annual one and the 14-month or Chandler frequency. From this point of 
view then, the determination of an appropriate interval for the computa- 
tion of the state is rather simple. The complications arise from the 
consideration of the capabilities that the satellite system can demon- 
strate. It is not only important that this system is accurate enough so 
that the signal can be separated from the noise; we must also be able to 
collect a sufficiently large number of observations over the estimation 
interval in order to be able to produce reliable results. The laser*sys- 
tems, being weather dependent, will have a disadvantage in that respect 
since observations lost on one day cannot be made up with additional ob- 
servations on the next day. That implies that a large number of optimal- 
ly selected stations should be employed at all times. 
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Because of the fact that the station polyhedron defines the CIS, a 
question which must be examined very carefully is the effect of missing 
stations' observations in the estimates. From the theoretical point of 
view, deleting or adding stations (or changing their location or coordi- 
nates) results in a new CIS, and one should therefore be alert to this 
for the proper interpretation of discontinuities or irregular behavior 
in the estimates. In connection with the reference system stability, we 
must also examine that of the satellite CIS. It is true that Lageos, 
which is the favored target for geodynamics research, has a very stable 
orbit compared to other geodetic satellites, but nevertheless we still 
do not know of an orbital model which would result in a centimeter level 
accurate orbit over a period of a few years. For obvious reasons it is 
not very desirable to have even a quasi-inertial system whose definition 
changes so much in so little a time. The general consensus on this issue 
is to use the satellite system as an interpolatory one and periodically 
calibrate it with respect to one of the systems that exhibits long-term 
stability such as Lunar Laser Ranging (LLR) and Very Long Baseline Inter- 
ferometry (VLBI). The application of the SRD mode, however, can consid- 
erably increase the stability of the satellite system due to the fact 
that the biases in the observations can be greatly reduced with proper 
scheduling. That not only will increase the quality of the end product, 
but it will also reduce the need for a frequent calibration procedure, 
which can be a nuisance and a source for further error. 

Returning now to the data distribution question, we should point 
out that by the discretization process we have in essence approximated 
the nonlinear functions which represent the coordinates of the pole with 
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a step function. It is obviously desired that the values of the step 
function at each interval be unbiased estimates of the average value of 
the true function over that interval. This can hardly be achieved if 
the collected data are not uniformly distributed over the entire inter- 
val as Fig. 24 illustrates. Even with a large number of stations in- 
volved, we still cannot guarantee the uniformity of the data at all times 
because we have no control over the weather or system failures. We sug- 
gest here that in processing data for polar motion determination, a dif- 
ferent set of coordinates of the pole should be estimated when the data 
set density changes abruptly, that is, each batch of observations which 
span a time interval no larger than that determined by the resolution of 
the systembeused to estimate one set of x,y coordinates. Furthermore, 
and perhaps more important, the reference epoch for these estimates 
should be computed on the basis of the data distribution rather than the 
middle epoch of the interval, as is currently assumed. In Fig, 25 we 
have plotted an assumed distribution of data for the problem of Fig. 24. 
With the current practice, the estimates will refer to epochs Mi, M 2 , 
etc. With the proposed new approach, the corresponding epociis will be 
El, E 2 , etc. When these epochs are used in plotting the variation curve, 
it is obvious that the estimated curve will be much closer to the true 
curve than the dashed line of Fig. 24 which is based on the Mi, etc. 
epoch labeling. The fact that the new estimates are not equally spaced 
should not alarm the standard users of this information, for a smooth 
curve can be fit to these points, and thereby one can obtain estimates at 
regular intervals. In fact, the astrometrically determined coordinates of 
the pole are always smoothed by means of Vondrak's method [1977] and 
the published results are indeed the output of this filtering process. 
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desired step-function approximation 

estimated curve 

— estimated step-function approximation 
observation intervals 

Fig. 24 The effect of nonuniform observing schedules on the 
estimation of the coordinates of the pole. 
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XvX.l Shaded area shows the quantity of collected 
observations over the indicated time interval 

Fig. 25 Unbalanced data distribution effects in the 
estimation of time-dependent averages. 
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5. NUMERICAL EXPERIMENTS AND RESULTS 
5.1 Simulation Studies 

A number of simulation studies were performed in order to substan- 
tiate our claims about the estimation of geodetic parameters from SRD ob- 
servations as opposed to the classical range observations. Simulations 
provide a lower bound on the accuracy of the results expected from the 
analysis of real data. In this case, however, this is of lesser concern 
to us because the main purpose of these simulations is to show the 
relative performance between the two approaches on the basis of identical 
data. 

In order though that the results of these simulations reflect real- 
ity too, we have used in most cases orbital models that contain all of 
the major perturbation sources and station configurations that either 
exist or whose existence in the near future is nearly certain. What has 
not been accounted for in the generation of the simulated observations is 
the weather. It is common practice in simulation studies to adopt a 
certain percentage p (usually 50%), which is used as a weather factor, 
i.e.jifN is the total number of observations which are possible, then pN 
of them are deleted to account for poor weather at the observing sites. 
Since we are interested in a relative comparison of the results, this is 
not an important issue. The treatment of the weather problem in the 
above fashion raises several questions, the most important of which is 
how does one apply the weather factor. It certainly makes little sense 
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to delete the first or last pN observations of the campaign since that 
only shortens its duration, and deleting pN observations at random is 
not very realistic either since it is equivalent to reducing the sampling 
rate by an appropriate factor. Modeling the weather realistically at 
each station, based on past weather records would certainly seem to be 
the most logical solution. This, however, would result in unduly in- 
creasing the complexity of the simulation process and with no major gains 
in the relative comparison of the results. 

In connection with the distribution of the data issue, a related 
factor is the observational capabilities of the stations. Although cur- 
rently only a few stations can observe during the day, most of the sta- 
tions are undergoing upgrading to the third generation of laser instru- 
mentation which will make it possible to observe at all times. For this 
reason, we have not discriminated between day and night observations in 
the simulations. 

In the following we discuss the three most important simulation 
studies performed. The first one involves nine existing stations. Its 
purpose was to find out what amount of SRD data could be collected from 
these nine stations over the observational period of ten days in the sec- 
ond half of August, 1980, and compare it with what was actually collected 
over that period of time. The purpose of the second simulation was to 
investigate the merits of using the proposed method in analyzing data 
collected from a hypothetical network of 17 stations which is likely to 
be realizable by 1983 as proposed in [ CSTG Bulletin , June 9, 1982]. The 
goal of this network is the establishment and maintenance of a Conven- 
tional Terrestrial Reference System by means of modern observational 
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techniques such as satellite laser ranging (SLR) and very long baseline 
interferometry (VLBI). Our study focuses on the problem of determining 
the baseline lengths between the observing stations of this network. In 
the third simulation, four nearly optimal baselines of the previous net- 
work have been singled out in order to study the capability of this sub- 
network in determining the coordinates of the pole. 

In all three simulation studies all possible simultaneous events 
between pairs of stations are generated. Since the SRD data set which 
consists of only simultaneous events is used to create the simple range 
data set, the number of observations in the latter is twice that of the 
former. In addition to this, because SRD observations will be in practice 
obtained from independent range observations, assuming that the two 
ranges in a pair have equal noise variance o^, the resulting SRD will 
have a noise variance equal to 2a^. This fact is also considered in the 
simulations. 

Each of the simulated pair of data sets (one for ranges and one for 
SRD's) is then used in estimating the parameters of interest, whether 
they be baseline lengths or the coordinates of the pole. Being simula- 
tion studies, if everything is left as is, there will be no difference in 
the recovered results larger than the input noise level. Our interest 
though is to compare the two methods when there are unmodeled orbital 
biases in the problem which are not being accounted for in the solution 
parameters. To achieve this, we apply a bias in the reference orbit which 
is a common input and identical for both the range and the SRD adjustment. 
The method which is least affected by this bias is obviously the one which 
recovers the parameters of interest closer to their "true" values. 
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5.1.1 Simulations for an existing SLR network. 

The fact that the simultaneity constraint reduces significantly 
the amount of usable ranging data in the SRD mode has been mentioned al- 
ready. In anticipation of using some existing laser ranging data for 
testing the proposed method, we performed an initial simulation study 
using the same stations from which the real data had been obtained and 
covering a period of time during which these stations seemed to have per- 
formed extremely well in acquiring large amounts of data. 

Of all stations observing during August of 1980, eight NASA sta- 
tions and two SAO stations seemed to be the only ones with significant 
amounts of data. After a preliminary inspection of the data distribution 
in time, it was apparent that nine of these stations had coobserved La- 
geos passes with significant overlapping intervals. Table 14 gives a 
summary of the available real data from all ten stations and the 


Table 14 Lageos Data Selection Summary 


Station Available Selected 

No. Observations Observations 


7090 

73390 

19042 

7943 

6068 

1934 

7092 

4167 

607 

7096 

4492 

786 

7120 

18245 

2749 

7907 

3040 

207 

7063 

5322 

1698 

71 15 

7163 

1047 

7091 

7322 

519 

71 14 

4130 

524 

Totals 

133921 

29113 


145 





corresponding amounts of data collected simultaneously byany two stations 
of the network. Comparing the totals in these two columns, it is obvious 
that less than 22% of the available data can be used in the SRD mode. 

From a detailed tabulation of the selected Lageos data (per day/ 
station), it was found that the greatest concentration is in the second 
half of the month for almost all stations. When the data collected in 
overlapping time intervals are selected, it turns out that only 38 
Lageos passes had been coobserved. A listing of the data per pass for 
each station in a station pair is given in Table 15. From this tabula- 
tion it becomes very clear that the two stations on the Australian conti- 
nent, Yaragadee (7090) and Orroral (7943), which dominate the complete 
data set (cf. Table 14), are also the ones with the most data in the 
selected (simultaneous) observations data set. The implications of these 
facts will become more obvious at a later stage when we will discuss the 
analysis of real data. From the initial 133921 available ranges, only 
2431 SRD events could be generated f-r the following seven station pairs: 
7090 - 7943 7115 - 7114 7114 - 7063 

7120 - 7115 7092 - 7120 

7943 - 7096 7092 - 7943 

Following the general guidelines of our simulation process, we have 
generated a data set of SRD observations which are collected from a net- 
work including the above eight stations and an additional one, the SAO 
station 7907, located at Arequipa, Peru. The generated data span a ten- 
day period, August 16 - 25, 1980, and we have used a sampling rate of one 
range every 30 s for all stations. A total of 26253 SRD events are possi- 
ble for the selected station pairs. The groundtracks of the observed 
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Table 15 Overlapping Data Distribution for Eight Station Pairs 


Sta t Ion 
No. 

Sa te 1 1 i te 
Paaa 

1 

1 

1 7090 

1 

V 1 

7943 

7096 7092 

7120 

7115 

71 14 

7063 

1 

1 664 







2 

1 



669 

647 



3 

I 351 

41 






4 

1 701 

76 






3 

1 



106 

30 



6 

1 836 

38 






7 

1 343 

43 






a 

1 431 

31 






9 

1 

13 

36 





10 

1 149 

32 






1 1 

1 117 

22 






12 

1 289 

34 






13 

1 220 

28 






14 

1 

6 

38 





13 

1 80 

a 






16 

1 640 

00 






17 

1 313 

12 






13 

1 313 

61 






19 

1 




935 

181 


20 

1 



66 

12 



21 

1 



286 

29 



22 

1 131 

10 






23 

1 100 

9 






24 

1 30 

6 






23 

1 




37 

40 


26 

1 




32 

26 


27 

1 



60 

37 



23 

1 226 

32 






29 

1 632 

88 






30 

1 176 

18 






31 

1 


6 

32 




32 

1 


1 12 

1 15 




33 

1 

10 

1 18 





34 

1 1217 

131 






33 

1 149 

46 






36 

1 876 

113 






37 

1 





228 

590 

33 

1 33 

7 
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Lageos passes and the baselines defined by the selected station pairs 
are shown In Fig. 26. We refer to this data set as the AG80 data, since 
the observational period falls In the Interval during which this campaign 
occurred. 

Even If we halve the number of possible events to account for the 
fact that not all of the stations had day-night observational capabili- 
ties, and further take half of the remaining observations to account for 
poor weather, equipment breakdowns, etc., we are still left with over 
6500 events. Considering that these events span only one-third of the 
entire month of August and, furthermore, that our sampling rate of one 
observation every 30 s Is nearly ten times lower than the observational 
rate of most stations, one should realize that the actual number of possi- 
ble SRD events for the entire month should be well above 50,000 with 
ample allowance for all factors. It therefore seems that either extreme- 
ly adverse conditions in the worst possible combinations dominated the 
performance of these stations over that period or there was a lack of 
effort in obtaining all possible observations. 

The data generated for the A680 data set were subsequently used in 
two recovery adjustments, one for polar motion components and one for 
baseline lengths. Both adjustments were done in the SRD as well as in 
the ranging mode. To identify which of the two modes is affected more 
severely by orbital errors, we perturbed the reference orbit in two dif- 
ferent ways. 

Initially we used only a random error applied to each coordinate of 
the satellite; subsequently, though, we augmented that with a linear 
trend time-dependent component. The adjustments are purely geometrical , 
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Fig. 26 Lageos groundtracks for the AG80 simulated data set. 


i.e., the reference orbit Is assumed perfectly known and no adjustment 
for It is allowed. The differences therefore between the simulated and 
recovered polar motion coinponents and baselines ("recovery errors") reflect 
the effect of the Induced orbital errors on the solved for parameters ot 
Interest. 

The random error used in the first simulation has a zprc mean and 
a 2.0 m standard deviation. The bias model of the second simulation can 
be analytically expressed as 

B. = a^ + b^(t - t^) + n. 1 = X.Y.Z (144) 

where the following (arbitrary) numerical values were used for the coef- 
ficients a^, b^ and the noise component n^; 

aj^ = -0.20 m bj^ = 0.01 m/day n^ - N(0.0, 0.05) 

By = 0.40 m by = -0.01 m/day 

a^ * -0.50 m b^ * 0.02 m/day (145) 

and t - t^ represents the time in days elapsed since the epoch t^ of the 

first day of the simulation. 

The recovery of the polar motion components was done for two differ- 
ent averaging interval scenarios. Initially, we broke down the ten-day 
mission into three intervals: a four day, a five day, and a one day, in 

this order. Subsequently, since the data distribution permitted it, we 
attempted a solution for ten daily averages. 

The baseline recovery simulation involved eight baselines between 
the same station network used for the polar motion simulation. Only the 
second orbit-biasing scenario was used in this simulation. 
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The results of the two simulations described so far have been sum- 
marized in the form of recovery errors in Table 16 for the polar motion 
components and Table 17 for the eight baselines. As can be seen from 
Table 16, the range adjustment results for the coordinates of the pole 
deteriorate by an order of magnitude when the orbit is perturbed by 
white noise (Case 1). In the case of the SRD adjustment though, the rms 
error increases by only 70% from its nominal value for the true orbit 
solution. When uhe errors in the orbit are of a systematic nature (equa- 
tion (144)) as in Case 2, the range adjustment results exhibit a fourfold 
increase in the rms recovery error compared to the SRD-based results. The 
results of the range adjustment can be significantly improved--al though 
still worse than the SRD results--by recovering the coordinates of the 
pole over shorter time intervals (Case 3). The less time the biases are 
allowed to accumulate, the less their effect on the estimated parameters. 
In this case the rms recovery error for the range adjustment dropped 
from 0'.'022 to less than 0'.'009 (Table 16, Case 3). 

The summary of the baseline recovery adjustments are given in Table 
17. Only the systematic error model was used in this simulation. The 
rms recovery errors were computed based oii the differences of the recov- 
ered baseline lengths from the a priori (modeled) baseline lengths. They 
are 1.9 cm and 3.7 cm for the SRD and the range adjustment results respec- 
tively. The contrast between these results is not as impressive as in the 
case of the coordinates of the pole, the reason being that we are dealing 
here (by choice) with extremely long baselines and in most cases with 
rather unfavorable geometry. In fact, from Fig. 26 one can gather that 
except for the baseline between stations 7943 and 7092, 
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Table 16 Polar Motion Component Recovery Error Summary 


Adjustment True Orbit Random Orbital 
Noise only Elrror 

Mode Case Case 1 


Systematic Orbital 
Elrror 

Case 2 Case 3 


Rnn,7e3 0.50 2.55 21.70 8.50 

SRD’s 0.40 0.69 6.10 6.20 


Note : Table values In ml 1 1 larcseconds . 


Table 17 Baseline Recovery Error Summary 


Baselines Adjustment 

between with 


J i n t i o U3 

Ran.'vcs 

SRD's 

7063 - 71 14 

-0.9 

-0.7 

7063 - 7907 

-3.2 

-3.0 

7090 - 79 43 

-i.a 

1 . 1 

7092 - 7120 

-2. 1 

-2.5 

7092 - 7943 

-3.3 

2. 1 

7096 - 79 -iS 

-3.2 

2..1 

7i iO - 7120 

0. 1 

0.5 

7120 - 79 -:3 

-0.9 

-1.3 

recovery error (cm) 

3.7 

1.9 


152 





in all other cases the geometry of satellite groundtracks - baseline 
direction is of the type to be avoided. This strong dependence of the 
SRD mode on the geometry of the problem has already been pointed out in 
Section 4.5.2, and these results merely confirm it. 

5.1.2 Simulations for a proposed SLR network. 

5.1.2. 1 Baseline Recovery . Based on an optimal global laser sta- 
tion distribution (likely to be realizable by mid-1983) proposed at a 
recent meeting of the COTES study group [CSTG Bulletin , June 9, 1982], a 
simulation study for baseline recovery was performed. Except for the 
fact that different stations (17 total) are involved, this simulation was 
similar to the one previously reported for the AG80 data set. The station 
locations and the data distribution are given in Tables 18 and 19. Base- 
line estimates and their statistics were computed for both the range and 
the SRD adjustments. In order to assess the effect of orbital biases on 
the baseline recovery, the orbit used in the adjustment (range and SRD) 
was again biased. In this simulation, however, the bias was applied in 
a slightly different manner from what was done in previous simulations, 
by applying it in terms of a radial, an along-track, and an across-track 
error as follows: 

radial bias 2.00 m 

along-track bias 0.60 m 

across track bias -1.20 m 

Two different adjustments were performed. In the first case the 
coordinates of all stations were obtained in a simultaneous adjustment 
based on the data collected from all station pairs. On the basis of this 
solution the baselines between all possible station combinations were 
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Table 18 Coordinates for the Stations Used in the 
Simulations 


Station ^ X (m) Y (m) Z (■) 


FTDAVS 7086 
WETTZE 7914 
HAWAII 7120 
STALAS 7063 
WESTFO 7091 
ONSALA 7095 
UERSTII 7911 
GRAZ 7999 
JAPAH 7935 
RICHIIO 7069 
QUIMCY 7051 
YARACA 7090 
CIIILBO 7901 
ORRORA 7943 
DIONYS 7940 
AREOUI 7907 
GRASSE 7942 


-1324510.442 
4074613.305 
-5464096.683 
1130304.818 
1492212.742 
3392750.872 
4022035 . 768 
4130031.490 
-4121637.800 
961333.601 
-2516274.896 
-2389125.331 
3844341.319 
-3912963.794 
4728637.251 
1941330. 1 15 
4530739.238 


-5332139.932 
931963.678 
-2402363. 153 
-4831721.449 
-4458121.791 
783278.257 
000000.000 
1106638.602 
3220176.370 
-5674186.968 
-4198843.469 
3042839 . 038 
-134247.357 
2259151.834 
1910493.462 
-5802024. 122 
639567.303 


3231791.056 
4801492.271 
2240358.273 
3993759 . 624 
4296005 . 489 
5325906 . 607 
4933350.635 
4716882.075 
363787 1 . 320 
2740519.741 
4075134.589 
-3078750.728 
5070549 . 690 
-4488060.975 
3817397.791 
-1796312. 9B6 
4408096.973 


Table 19 Distribution of Ranges and SRD's for 
Each Baseline 


No. 

Base line 
End Stations 

Lenf^th 
( m) 

Ran^e 

Obs. 

SRD 
Obs . 

1 

7901 

= > 

7914 

1123487.006 

7202 

3601 

2 

7095 

= > 

7940 

2308833.694 

5976 

2980 

3 

7942 

= > 

7999 

700368. 121 

7460 

3730 

4 

7095 

= > 

7942 

1484391.097 

6954 

3477 

3 

7091 

= > 

7093 

5669637.481 

2594 

1297 

6 

7063 

= > 

7911 

5708839 . 39 1 

2442 

1221 

7 

7069 

= > 

7942 

7431634.061 

446 

223 

8 

7911 

= > 

7940 

2322728.588 

5884 

2942 

9 

7901 

= > 

7942 

1239620.646 

7062 

3531 

10 

7942 

= > 

7914 

683352.290 

7468 

3734 

1 1 

7911 

= > 

7095 

1073641.514 

7450 

3725 

12 

7942 

= > 

7940 

1412734.344 

6496 

3248 

13 

7093 

= > 

7999 

1009482.790 

7478 

3739 

14 

7999 

= > 

7940 

1346693.332 

6620 

3310 

13 

7093 

= > 

7914 

872957. 1 16 

7612 

3806 

16 

7091 

= > 

7069 

2044497 . 683 

4496 

2248 

17 

7063 

= > 

7907 

5926366.472 

1024 

512 

13 

7036 

-> 

7907 

6014011.633 

994 

497 

19 

7069 

= > 

7907 

4643137.992 

1716 

858 

20 

7069 

= > 

7036 

2363121.040 

4112 

2056 

21 

7063 

= > 

7031 

3701936.397 

3944 

1972 

oo 

7031 

= > 

7036 

1 043222. 236 

5296 

2648 

23 

7120 

= > 

7051 

3909408. 132 

2950 

1475 

24 

7120 

= > 

7026 

5167466.031 

1830 

925 

23 

70E5 

= > 

7063 

2618612.747 

4666 

2333 

26 

7091 

= > 

7036 

3135343.207 

4348 

2174 

27 

7120 

= > 

7935 

5947116.046 

1342 

671 

23 

7933 

= > 

7090 

7171939.095 

366 

183 

29 

7090 

= > 

7943 

420CO79 . 994 

3226 

1613 

30 

7935 

= > 

7031 

7603305.998 

412 

206 
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% 

obtained along with their fohma^ccuracies and differences with respect 
to their "true" values. Part of the baseline results are summarized in 
Table 20. Only the cases in which a station pair has coobserved have 
been listed. In all cases but two, one listed (7090-7943) and one not 
(7901-7911), the baseline lengths have been overestimated although the 
errors in the SRD case are about an order of magnitude smaller than the 
ones for the range adjustment. Since the positive radial bias results in 
an "expansion" of the network of satellite positions, this should come as 
no surprise. The stations have a global distribution and because the 
observations from all stations are adjusted simultaneously, their posi- 
tions become interdependent and the aforementioned expansion affects all 
of them similarly. Fig. 27 shows the results of the two adjustments for 
all possible baselines. 

This first solution prompted us to test the recovery of baselines 
from independent adjustments. In this second case the data collected 
from each pair of stations are adjusted independently and the estimated 
baselines are only the ones defined by coobserving station pairs. The 
results of this second type of solution are shown in Table 21. What is 
obvious again is that the SRD results for the baseline lengths are again 
superior to the range adjustment results. 

The mvtst interesting observation, though, in this solution is that 
on the basis of the same data the range adjustment now underestimates the 
baselines and the recovery errors are all negative. The reason behind 
this is the one-sided data distribution in this instance, as opposed to 
the global distribution that existed in the network adjustment case. 

From Fig, 26 it is obvious that we are dealing with extremely long 
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Table 20 Baselines Recovered from Network Adjustments 
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RECOVERY ERRORS 




BASELINE LENGTH - (KM) 


Fig. 27 Network solution results, recovery errors versus 
baseline length. 


Table 21 Baselines Recovered from Independent Adjustments 
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baselines in which case the simultaneous events are confined in the area 
in between the end stations. The loss of the uniform data distribution 
around the stations results in a one-sided biasing of the station posi- 
tion towards the opposite station or, better, towards the "barycenter" 
of the observed satellite events. Since this, as we discussed above, 
lies in between the two stations, thence the "shrinkage" of the estimated 
baselines. 

For the SRD results there seems to be no bias preference, and 
those errors are rather randomly distributed and in almost all cases at 
the centimeter level. The three baselines for which the range adjustment 
has given better results than the SRD all have lengths in excess of 7000 
km and very few observations. As it has been previously reported, the 
SRD mode is much more geometry dependent than the range mode, and as the 
results of Table 21 show it admits of its limitations very eagerly (note 
the formal accuracies on those baselines!). Unlike the SRD mode, the 
formal accuracies for the range mode give no hint whatsoever as to the 
real accuracy of the results. Even though the recovery errors are of the 
order of a few decimeters in all cases, the reported a's are hardly ever 
higher than 2 cm! A pictorial presentation of the recovery errors for 
this solution are shown in Fig. 28. 

On the basis of these simulations one can conclude that the SRD 
mode will in all likelihood provide more meaningful results in the pres- 
ence of unmodeled orbital biases of the type considered herein than the 
range mode, and it will also give more reliable accuracy estimates for 
those results. Comparing the batch (global) solution to that of individ- 
ual adjustments, the latter seems to be by far a better approach in the 
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Fig. 28 Independent solution results, recovery errors versus baseline 
length. 
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Table 22 Summary Statistics for Baseline Recovery Errors 


Adjastment Node 

; 

Network 

Independent 

Stat lat Ic 

Obaervab le 



Hean 

Recovery 

Ranges 

89 

-55 

Error (cm) 

SRD* a 

18 

4 

Mean 
Ra t lo 

Rani^ea 

S2 

-17 

Error/Stgma 

SRD’ a 

17 

0.2 


case of SRD observations, although the opposite is true for the range 
observations. Compare, for instance, the level of recovery errors be- 
tween Tables 20 and 21. This is also documented by the average values 
of the recovery errors displayed in Table 22 for both the ranges and the 
SRD's as obtained from the network as well as the independent solutions. 

5. 1.2.2 Polar Motion Parameter Recovery . Since the global rota- 
tions of the CTS polyhedron will be monitored by a subset of the defining 
stations [ CSTG Bulletin , June 9, 1982], we selected four baselines out 
of the 136 possible between the proposed 17 stations [ibid.] to investi- 
gate the performance of sucn a subnetwork in estimating the coordinates 
of the pole. 

The selected baselines were chosen in such a way that they conform 
as nearly as possible with the optimality criteria established in the 
previous chapter of this investigation with respect to polar motion esti- 
mation from SRD observations. The locations of the eight stations defin- 
ing these baselines are given in Fig. 29. A set of SRD observations was 


161 





ORIOtNAl' PACw ’S 
OF POOR QUALITY 



ring. 




ORIGINAL PAGE IS 
OF POOR QUALITY 

generated, spanning a seven-day period, August 14 through 20, 1980. The 
sampling rate in this simulation was one observation per minute, and to 
each SRD observation was added a white noise component with a standard 
error of 14 cm, the equivalent of 10 cm in each range. The coordinates 
of the pole were specified as two-day averages; their "true" numerical 
values are given in Table 23. The satellite orbit was biased in the same 
manner as for the previous simulation, this time by the following errors: 

radial bias 1.00 m 

along-track bias 0.06 m 

across-track bias -0.12 m 

Table 23 Polar Motion Component Values Used 
in the Simulations 


Tlinc 

Interva 1 

True ' 

Va lues 

YYIEiDD 

riTIilDD 

X 

y 

S003K 

- 800316 

-0"020 

0"310 

CJOQK' 

- 300313 

-0"019 

0"31 1 

300313 

- 300320 

-0"013 

0''312 


The results of the two adjustments for the three two-day averages 
of the coordinates of the pole are shown in Table 24. The rms errors for 
the range adjustment are nearly an order of magnitude higher than those 
for the SRD adjustment. Table 25 lists the statistics of the estimated 
coordinates of the pole for both adjustments. Despite the comments made 
above as to the quality of each solution, the formal statistics give no 
hint at all about it, and it should once again be pointed out that they 
are completely unreliable in the presence of unmodeled errors. 
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Table 24 Polar Motion Component Recovery Results, Complete 
Data Set Solution 


simulated Pole Coordinates 


Recovered Pole Coordinates 


Time 

iuterva 1 

True Values 

Raume Solution 

SRD Solution 

YYIillDD 

Y^TIIIDD 

X 

y 

X 

y 

X 

y 

B00314 

- 300316 

- 0"020 

0*310 

-0*008 

0*316 

-0*019 

0*308 

000016 

- 800.11 a 

-0*019 

0*311 

-0*004 

0*317 

-0*019 

0*311 

eoosiG 

- 300320 

-0*018 

0*312 

-0*006 

0*318 

-0*017 

0*31 1 

rino recovery error 


0*013 

0*006 

0*001 

0*001 


Table 25 Statistics for the Recovered Polar Motion Components 
Obtained from the Complete Data Set 



P. fi n ,7 

e So 

1 u t Ion 

S R 

D Sol 

u t 1 o 

n 

Step 

Ctandar J 

De V 1 a 1 1 o 113 

Corre lu t ion 

Standard 

Devla t Iona 

Corre 

la t Ion 

( 1) 

0’:r( 1) 

ffyl 1) 

x( 1) <-> y( 1) 

Cxi 1) 

Cyl 1) 

x( 1) < 

-> y( 1) 


( !!!a8 ) 

( nna) 

( mas) 

( mas ) 


1 

0.0S9 

0.083 

-0.023 

0.314 

0.459 

-0. 

148 

2 

0.039 

0.0S7 

-0.046 

0.31 1 

0.459 

-0. 

142 

o 

0.090 

0.037 

0.003 

0.311 

0.459 

-0. 

145 


It must be obvious by now that an order of magnitude Improvement in 
the accuracy of the coordinates of the pole can be achieved by analyzing 
the exact same range observations in the SRD mode. One should also con- 
sider that this particular simulation is based on feasible station loca- 
tions which are in no way the optimal network for polar motion determina- 
tion, and additionally no effort was made to single out the optimal sat- 
ellite passes for each specific baseline. Yet, the results are rather 
promising in view of the generally accepted requirement of O'.'002-accurate 
pole positions over two-day intervals. 

In an effort to account for the effect of loss of data we have 
readjusted the simulated observations only this time we restricted the 
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admissible observations to those with an elevation of 40° above the 
station's horizon. The recovered coordinates of the pole and their sta- 
tistics are given in Tables 26 and 27 respectively. The SRO adjustment 
still recovers the parameters better than the range adjustment. What is 
more interesting here, however, is the indication of a downward trend in 
the correlations between the recovered values for the SRD case and the 
opposite for the range case. When we used the more restricted data set 
in the second simulation, we in effect forced our observations to lie 
within and around the two baseline stations. As we have already seen in 
Section 4.5.3, this is the optimal region for selecting observations 
which will result in the highest sensitivity with respect to one of the 
parameters and the lowest with respect to the other. Although the 


Table 26 Polar Motion Component Recovery Results, Restricted 
Data Set Solution 


Siuiilated Pole Coordinates 
Tina Interval True Values 


VYIIiIDD 

7L 

y 

800314 - SOOSlu 

-0"020 

0"310 

C0C31O - aoosia 

-0“019 

0“31 1 

300313 - 300320 

-0-013 

0"312 

i’uiJ recovery error 



Recovered Pole Coordinates 
Raupre Solution SRD Solution 


K 

y 

X 

y 

0"001 

0"31 1 

-0"021 

0*308 

0"04)5 

0"305 

-0"O19 

0*310 

0“0I7 

0*303 

-0"0I9 

0*309 

0 "028 

0"007 

0*000 

C“002 


Table 27 Statistics for the Polar Motion Components As Obtained 
from the Resti-lcted Data Set 



R a u ,7 

c ? o 

1 11 t ion 

S R 

D Sol 

u ! 1 o n 

S' ep 

Siaudard 

9e via 1 io*ic 

Corre la t ion 

Standard 

De V 1 a 1 1 o ns 

Corre la t ion 

V i) 

0 -:< i) 

Gy( 1) 

3( 1) <-> y< i) 

Gx( 1) 

Cy( 1) 

x( 1) <-> y( 1) 


( iiias ) 

( iiias ) 


( mas ) 

( mas ) 


1 

0.43-1 

0.597 

o.ooa 

1.203 

1 . 161 

-0. 103 


0.493 

0.613 

-0.016 

1 .219 

1.112 

-0. 104 


0.493 

O.oOo 

0.006 

1.246 

0.955 

-0.044 
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original correlation of -15% is not really significant, it is only fair 
to point out that in the restricted second simulation it drops to -10% 
or less. At the same time the range adjustment results show in general a 
slight Increase in the correlation. 

Finally, we compare the formal accuracy estimates between the two 
solutions. The first simulation is based on 7407 SRD's, while the second 
on only 669, a ratio of about 11:1. Since in both cases the noise level 
is the samej one would expect that the increase in the formal accuracies 
between the two adjustment should follow the ^ law, n being the ratio of 
the observations between the two solutions, or in our case we would ex- 
pect roughly a threefold increase in the o's of the recovered parameters. 
Any deviation from this ratio on the higher side indicates that the 
poorer data distribution affects the solution, while on the lower side it 
indicates that the new geometry is superior to the former. This is ex- 
actly what happens liere since the a-ratic for the ranges is nearly 5.6 
(compared to the expected /I 1.07 s 3.3), while for the SRD's it is only 
2.5. The invoked data selection has not only compensated for the loss 
of data, but in fact it has improved the sensitivity of the system with 
respect to the estimated parameters. 
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5.2 Experiments with Real Data 

5.2.1 Preliminary adjustment for data editing. 

Since the purpose of this investigation is the introduction and 
study of a new method for the analysis of ranging observations to satel- 
lites and in particular to Lageos, knowing the quality of the data which 
we planned to use for testing the method was one of our primary concerns. 
As we have already discussed, the distributing agencies such as NASA 
and SAO do edit the raw data and delete most of the spurious observations 
This editing, however, is done on a pass-by-pass basis and not for 
the whole aggregate of the available data. To do so one has to "fit" an 
orbit to the data as a whole and to compare the discrepancies of each 
observation from that orbit. This is the primary reason for which we 
considered an adjustment for the complete data set as a necessity. 

Secondarily, though, such an analysis of the data would also pro- 
vide an indication of how well and to which level of accuracy our orbi- 

tal model fit the data. The lower the rms residual of the observations, 
the better the orbi cal model deployed. The qualitative and quantitative 
characteristics of this model are of great interest in this investigation 
We will later use this model to calculate the reference orbit with re- 
spect to which the SRD observations are adjusted. 

Last but not least, since most of the data have already been ana- 
lyzed by other agencies, the results of our own software (GEOSPP— GEO - 

detic Satellite Positioning program) could provide a check on the code 
and give us some confidence in the program. This latter is always a 
major problem since, as someone put it, "Every nontrivial program has at 
least one error." 
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An initial inspection of the August, 1980, Lageos data set (Sec- 
tion 5.1.1) had already indicated that the highest concentration of si- 
multaneously observed Lageos passes fell in the second half of the month 
of August. For that reason our effort focused on the subset of data 
spanning the period from August 14 through 31, 1980. It turns out that 
the last coobserved pass occurs on August 28. Data from the ten most 
active stations during that period were selected in order to subsequent- 
ly be used for SRD event generation. A total of 24240 range observations 
were selected in such a way that all stations but one (7907) have nearly 
the same amount of data. These data were analyzed with our computer 
program GEOSPP in a preliminary adjustment. 

The theory behind the orbital model used in GEOSPP has been de- 
scribed in Section 2.2; the numerical values for some of the constants 
used in the program are given in Table 28. We have used the geopotential 
coefficients of the preliminary model PGS-Ll [Lerch and Klosko, 1981] up 
to degree and order twelve, since as it is reported in [ibid.] the per- 
turbations of higher harmonics for such a short period of time (14 days) 
are nearly equal with the errors caused by the uncertainties in the 


Table 28 Numerical Values of Constants Used by GEOSPP 


Semi-axis major 

Inverse flattening...... 

Gravitational constant.. 

Rotational rate 

Speed of light 

Astronomical Unit (i AU) 
Solar pressure at 1 AU. . 


6378144.11 

( m) 

298.255 

— 

398600.4125x10^ 

(mVo2) 

0.000072921158547 

( rad/a ) 

299792458.0 

( la/a) 

149547870950.0 

( ra) 

4.62576x10* 

(H/m2) 
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coefficients themselves. An initial state vector for Lageos was pro- 
vided by Dr. P. Dunn, although in a reference frame different from that 
used in GEOSPP. The state vector that we used was obtained from the 
former by applying to the former the nutation, precession and equinox 


Table 29 Numerical Orbit Integration Information for GEOSPP 


VARIABLE ORDER- VARIABLE STEPSIZE RUMERICAL INTECRATIOR : 


140.0 SEC 
0.0 SEC 

600.0 SEC 

1 . 0D-07 
1 . 0D-04 
21 


PERTURBATION MODEL : 


NOMINAL STEPSIZE : 

MINIMUM STEPSIZE s 

MAXIMUM STEPSIZE : 

RELATUT ACCURACY 

EQUATIONS OF MOTION : 

VARIATIONAL EQUATIONS...: 
INT.rCSSACE OUTPUT UNIT : 


GEOPOTFJITIAL : (12,12) 

FOR VARIATIONAL EQUATIONS : ( 4, 4) 

MOON : YES 

SUN : YES 

UNIIODELED ACCELERATIONS : YES 

ALONG TRi\CK = -0.3480D-11 (M/S2) 

CROSS TRACK = 0.0 (M/S2) 

RADIAL = 0.0 (M/S2) 

SOL.AR RADIATION PRESSURE YES 

SATELLITE AREA = 0.2827 (M2) 

SATELLITE MASS = 406.9650 (KG) 

SATELLITE REFLECTIVITY CR= 1.1729 

SOLID EARTH TIDES YES 

LOVE NO. !<2= 0.2740 

PHASE ANCLE £2=2.3300 


LOVE NUMBER FOR RADIAL EXPANSION H2 = 0.600 

SUIDA NUMBER FOR UCRIZONTAL SHEAR L2 = 0.075 


INPUT INFORILATION FOR ARC : 7603901.01 


YYMMDD HHMMSS.SSS 

EPOCH OF ELEMENTS : 80 813 235930.028 

OBSERVATIONS START AT. .: 80 814 00 0.000 

OBSERVATiOI.S END AT : 80 829 0 0 0.000 

REFERENCE SYSTEM EPOCH.: 80 731 


INERTLAL CARTESIAN ELEMENTS AT THE EPOCH : 


X (ID 

-5398161.430 

XDOT (?I/S) 
-4720.0626632 


Y (M) 

-5961700. 182 

YBOT (M/S) 
-870.0206757 


Z (M) 

-9377042.063 

ZDOT (M/S) 
311 1.3650676 
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corrections that relate the two reference frames. The result along with 
the other information required for the numerical integration of the 
orbit is given in Table 29. The orbit was allowed to adjust freely in 
this solution. The a priori station locations are given in Table 30. 
Table 31 gives the summary of the residuals' statistics for the last 
iteration of the solution. We have already discussed the fact that the 


dynamic mode is an ill-conditioned problem due to the physical 


Table 30 A Priori CTS Station Coordinates 


Stat Ion 

No . 

X (n) 

Y (m) 

Z ( d) 

STALAS 

7063 

1130711.700 

-4831371.300 

3994088.700 

ML05O2 

7090 

-2389002.297 

5043333.488 

-3078528.544 

10.0702 

7091 

1492450.900 

-4457281.700 

4296817.000 

HLO802 

7092 

-6143448.500 

1364706.900 

1034164.800 

10.0602 

7096 

-6100049.584 

-996197.831 

-1568978.317 

F0.02 1 1 

7114 

-2410428. 196 

-4477802.221 

3838688.071 

MLO307 

71 15 

-2350B67.357 

-4655546.092 

3660999 . 228 

F0.0110 

7120 

-5466003.686 

-2404404.305 

2242228.593 

ARELAS 

7907 

1942786. 100 

-5804078.900 

-1796938.600 

ORRLAS 

7943 

-4447545.681 

2677137.812 

-3694997.951 


Table 31 Residual Summary for the Complete Lageos Range Data 
Set Adjustment by GEOSPP 


ADJOSTHEirr STATISTICS FOR ITERATIOIf : 2 


DECREES OF FREEDOM FOR THIS ADJUSTMENT : 24234 

PREVIOUS WEIGHTED SUM OF SOUARES OF THE RES ID.IALS/D. F : 0.1760 

CUIWENT WEIGHTED SUM OF SQUARES OF THE RES IDLALS/D. F : 0.1760 

IMI'ROVEMENT IN PERCENT (NEC. SIGN INDICATES DECREASE) : 0 

CONTRIBUTION FROM STATION PARAMETER CONSTRAINTS : 1745. D-04 

CONTRIBUTION FROM POLAR MOTION PARAMETER CONSTRAINTS : 5508. D-12 

CONTRIBUTION FROM ORBITiU, PARAMETER CONSTRAINTS : 0. 


PASS BY Pi\SS BREAKDOV/N OF ADJUSTMENT STATISTICS FOR ITERATION s 2 


PASS CONSTRAlIflS FROM : TOTAL NO. OF NUMBER OF NUMBER OF 

NO. STATIONS P.M.STEPS ORBIT CONSTRAINTS OBSERVATIONS PARAMETERS 
I 30 60 36 24240 42 


WEIGHTED SS. DEGREES OF 
OF RES I DUALS FREEDOM 
0.42666D+04 24234 


VARIANCE VITV(I) RMS OF THE MEAN OF 
COMPONENT #0BS.-6 RESIDUALS RESIDUALS 
0.1761 0.18 0.4195 -0.0009 
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characteristics of the earth (Section 4.3). This ill-conditioning can 
be alleviated in a range data solution if, for example, the longitude of 
one station is constrained. In this particular problem though we prefer- 
red an alternate solution which has been discussed previously in [Pavlis, 
1979], that is, the "ridge estimation." By applying a very small weight 
in all station coordinates, we overcome the numerical singularity of the 
normal equation matrix and at the same time we show no preference for 
any single station. The weight applied in this case corresponds to a 
variance of (50 m)^ in each station coordinate. 

As the residual summary in Table 31 shows, the data seem to be of 
rather good quality since they fit the orbit with an rms of 42 cm over- 
all. The station position estimates and their statistics are given in 
Table 32 and the adjusted satellite state vector and its statistics in 
Table 33. It should be no surprise that the estimated variances for the 
X and Y coordinates are so much higher than those for Z. The ridge esti- 
mator has overcome the numerical singularity, but it does not separate 
completely the two parameters which are associated with the stations' 
longitude, i.e., X and Y. This should be of no inmiediate concern, since 
the relative quantities such as the baseline lengths which are of more 
interest to us are not affected by this peculiarity. 

The estimates for the baseline lengths for all possible station 
pairs (45 total) are given in Table 34 along with their formal accuracies. 
All of these estimates show an internal consistency of 2 - 3 cm. Excep- 
tions are the baselines which involve station 7907, primarily because 
that station has almost ten times less data than the average station in 
the solution. 
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Table 32 Station Coordinates and Standard Deviations Estimated 
by 6E0SPP 


■ tatUa 

Ra. 

X (m) 

Y (■) 

Z (■> 

OX(B> 

CSY(b) 

OZ(b) 

0TALAB 

TM3 

1IM7M.44S 

-46JI37I.26I 

3994*66.938 

8.726 

1.341 

* 

*4* 

mj$9»a 

7«9« 

-22B4— a 

8646.;?4.aSI 

-3*76827.412 

8.979 

2.633 

* 

*38 

HL»7»a 

7MI 

14*2446.214 

-4487262.4*4 

4294617.648 

8.264 

1.77* 

* 

*41 

HL*M2 

7M2 

-4143446. 162 

13447*9. 1 18 

1*34143.844 

1.416 

7.264 

* 

*41 

HLM*2 

7M« 

-6I46649.460 

-994198.2*2 

-1846977.313 

1.161 

7.232 

* 

*41 

NL2211 

7114 

-34I443*. 112 

-44770*1.264 

3636466.977 

8.3*9 

2.636 

* 

*4* 

NL9M7 

7IIS 

-239666«.2I6 

-4488848. 147 

344*999.686 

8.819 

2.787 

* 

*37 

NL«II* 

712* 

-B466M4.676 

-24*44*2. *6* 

2243329.441 

2.681 

4.461 

* 

•44 

ARKLAS 

7*#7 

1942763.226 

-86*4*6*. 2*8 

-1794919.713 

4.661 

2.3*4 

* 

*91 

0RRLA8 

7*42 

-4447344.472 

2477 14*. 227 

-3494997.474 

3. 174 

8.273 

* 

*34 


Table 33 Initial State Vector for Lageos As Obtained by GEOSPP 
from the Complete Lageos Range Data Set Adjustment 


Reference Syaten : 


Mean of 800731 


POSITION 

X (m) 

Y (m) 

Z (m) 

Es t Ima te : 

Standard Deviation : 
rma position ( n) 

-9097881.421 

7.068 

-9961549.913 

6.044 

5.369 

-9377290.436 

0.042 

V E L 0 C I T Y 

X ( m/s ) 

Y (m/s) 

Z ( m/s ) 

Es t ima te 

Standard Deviation : 
rms velocity (m/s) : 

-4720.07387 

0.00103 

-87 1 . 25438 
0.00560 
0 . 00329 

3111. 28230 
0.00002 

Reference System : 

TOD 

800813 235930. 

028 

Pos i t ion ( m) : 

Ve loc 1 ty ( m/s ) : 

-5098243.471 

-4720.05663 

-5961677.201 

-870.98060 

-9377012.675 
31 1 1.38509 


Note ! TOD - "True of Date" 
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Table 34 


POOR 


PAGE IS 


Baseline Lengths and Standard Deviations Estimated 
by GEOSPP 


Base 1 Ine 
No. 

End 

Sta t Iona 

Apr lor 1 
Ea t Ima te 
( n) 

Adjun ted 
Value 
( n) 

Adjua ted 
ininua 
Apr lor 1 

Sienna 
( n) 

1 

7063 

« = > 

7090 

12645981.761 

12645950.847 

-0.914 

0.018 

2 

7063 

= = > 

7091 

602032. 143 

602032. 169 

0.026 

0.036 

3 

7063 

= »> 

7092 

10003296.513 

10003295.833 

_0.682 

0.023 

4 

7063 

= = > 

7096 

9896473.055 

989647 1 . 526 

- 1 . 528 

0.022 

3 

7063 

= = > 

7114 

3562138.713 

3562137.442 

- 1 . 272 

0.041 

6 

7063 

= a> 

7113 

3301893. 178 

3501891.797 

-1.381 

0.037 

7 

7063 

= = > 

7120 

7244020.742 

7244019.261 

-1.482 

0.028 

8 

7063 

= = > 

7907 

5928036.931 

5928019.003 

-17.948 

0.085 

9 

7063 

= = > 

7943 

12108339.034 

12108538.064 

-0.990 

0.018 

It) 

7090 

= »> 

7091 

12638160.062 

12638160.219 

0. 156 

0.024 

11 

7090 


7092 

6674009 . 770 

6674008.743 

- 1 . 027 

0.024 

12 

7090 

= = > 

7096 

7247520.432 

7247520.743 

0.311 

0.023 

13 

7090 

= = > 

7114 

11768618.014 

11768618.337 

0.323 

0.018 

14 

7090 

= = > 

7113 

11810628.856 

1 1810629.014 

0. 158 

0.016 

13 

7090 

ss> 

7120 

9656458.379 

9656438.910 

0.331 

0.021 

16 

7090 

- = > 

7907 

11730456. 119 

11750458.620 

2.500 

0.034 

17 

7090 

S5> 

7943 

3196328.733 

3196328.646 

-0.087 

0.021 

18 

7091 

= = > 

7092 

10141371.223 

10141371.602 

0.379 

0.031 

19 

7091 

= = > 

7096 

10199643. 124 

10199642.536 

-0.587 

0.025 

20 

7091 

3S> 

7114 

3929728.800 

3929728.019 

-0.782 

0.039 

21 

7091 

= = > 

7113 

3900398.443 

3900397.570 

-0 . 876 

0.034 

22 

7091 

= = > 

7120 

7540273.824 

7540273. 123 

-0.701 

0.029 

23 

7091 

= = > 

7907 

6257037.782 

6257020.271 

-17.511 

0.088 

24 

7091 

= = > 

7943 

12249596.212 

12249596.272 

0.059 

0.022 

23 

7092 

= »> 

7096 

3514556.686 

3514334.371 

-2.316 

0.027 

26 

7092 

= = > 

7114 

7479017.596 

7479018.461 

0.865 

0.027 

27 

7092 

= = > 

7118 

7384680.410 

7584681. 155 

0.745 

0.028 

28 

7092 

= = > 

7120 

4015538.430 

4015538.979 

0.550 

0.028 

29 

7092 

= = > 

7907 

11171113.715 

11171110.424 

-5.291 

0.041 

30 

7092 

= = > 

7943 

5192643.026 

3192640.982 

-2.044 

0.024 

31 

7096 

= = > 

7114 

7414696.951 

7414696.912 

-0.040 

0.023 

32 

7096 

= = > 

7113 

7402692.901 

7402692.731 

-0. 170 

0.023 

33 

7096 

= = > 

7120 

4112220.342 

4112220.461 

-0.081 

0.025 

34 

7096 

S = > 

7907 

9373094.032 

9373093.497 

-0.554 

0.044 

33 

7096 

= a> 

7943 

4534571.701 

4554572. 165 

0.464 

0.824 

36 

7114 

= s> 

7113 

258289 . 958 

258290. 167 

0.210 

0.036 

37 

7114 

= = > 

7120 

4022959.527 

4022959.505 

-0.022 

0.031 

38 

7114 

= = > 

7907 

7243602. 178 

7243388.024 

-14. 134 

0.076 

39 

7114 

= s> 

7943 

10587702.281 

10587702.701 

0.420 

0.018 

40 

7113 

s = > 

7120 

4096904. 174 

4096904. 146 

-0.027 

0.031 

41 

7113 

= = > 

7907 

7038726 . 637 

7038712.246 

-14.411 

0.074 

42 

7115 

= = > 

7943 

10595990. 172 

10595990.425 

0.253 

0.018 

43 

7120 

ss> 

7907 

9097407.601 

9097399 . 389 

-8.212 

0.052 

44 

7120 

= = > 

7943 

7880988.899 

7880989 . 300 

0.401 

0.022 

43 

7907 

= = > 

7943 

10787493.058 

10787496.733 

3.676 

0.041 


The station-by-station and pass-by-pass analysis of the residuals 
gives some more insight into the relative performance of the stations and 
the relative quality of their data. Tables 42 through 51 in Appendix D 
give these sun.maries for each of the ten stations. In general, the stan- 
dard deviation of the residuals in a pass is at the decimeter level, 
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although for some of the stations (e.g., 7063f 7907, and 7943) it is 
about three times higher than that or more. In addition to that, it 
seems that the observations from three of the stations on four particu- 
lar occasions include an unreasonably high bias. These are 
for station 7063: pass #4 -2.1 m 

(cf. Table 42), 

pass #8 -4.4 m 

for station 7114: pass #11 -4.5 m (cf. Table 47), 

and for station 7115: pass #4 -6.8 m (cf. Table 48). 

Overall, the data seem to be of consistent quality, except for the in- 
stances cited above, and the performance of the orbital model and the 
computer software were satisfactory. 

The data were subsequently examined to find the baselines with the 
most simultaneously observed passes. A computer program (OVERLAP) sup- 
plied to us by Mr, R. Kolenkiewicz of NASA/GSFC was used for this purpose. 
When we isolated the data falling in the overlap periods and examined 
their distribution by station and by time, it was realized that there 
were no aggregates of passes that spanned intervals of time long enough 
to detect polar motion with decimeter level observations. In addition to 
that, the number of observations per baseline was disappointingly small 
to attempt a baseline solution, except perhaps for the station pair 
7943-7090 which had 984 observations. For these reasons we concentrated 
on attempting a solution in the SRD mode with data collected from the 
aforementioned baseline only. 
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5.2.2 Estimation of the 7943-7090 baseline. 

The overlap range data from stations 7943 and 7090 were processed 
along the guidelines established in Chapter 3 for the generation of 
"quasi-observable" SRD's. The data distribution for the range observa- 
tions from each station is given in Tables 35 and 36. Out of these 
nearly 7500 range observations, only 975 simultaneous events could be 
generated. The corresponding range data set of course contains exactly 
twice as many observations, i.e., 1950. 


Table 35 Observation Summary for Range Data from Station 7943 


Station : 7943 PaBS«>a 'i'racked : 32 Obaerva t Iona total s 3418 

s = ssss3srssssssa3rsszzs::: = sssr:ssssssrBs33sarssrs s sss-s szsssbsbsxs ssbss s b 3 s 


Pass 
Wo . 

Bcalnnlngr Dr. te 
YYMMDD HHMMSS.S 

End Ing 

yymmdD 

Da te 

HHHltSS.S 

Paaa 

Dura t Ion 
(a) 

Oba. 

Dens 1 ty 
Laff 
( a ) 

1 

800814 

1191 0.2 

800BI4 

1 12015.2 

1 154.9 

56 

20.62 

2 

800813 

92743 . 0 

800813 

100630.3 

2323 . 3 

122 

19.06 

3 

800815 

124959.9 

800815 

132145.0 

1905. 1 

82 

23.23 

4 

800815 

1603 7.7 

800813 

1641 0.0 

2272.3 

98 

23. 19 

3 

800813 

1943 7.6 

800815 

200932.4 

1484.7 

83 

17.89 

6 

800816 

81544.9 

800816 

84152.7 

1567.8 

85 

18.44 

7 

800816 

114022.8 

800816 

1204 7.8 

1425.0 

90 

15.83 

8 

800816 

144522.7 

800816 

1531 7.7 

2745 . 0 

212 

12.95 

9 

800816 

181315.0 

890816 

190345.0 

2910.0 

159 

18.30 

10 

800817 

101632.8 

800317 

1049 7.7 

1934.9 

47 

41. 17 

11 

80081V 

134930.4 

800817 

140545.0 

1514.7 

54 

28.05 

12 

800817 

170113.0 

800817 

1728 0.0 

1605.0 

75 

21.40 

13 

809818 

90422.7 

BOOH 18 

93445.2 

1822.5 

139 

13. 1 1 

14 

800818 

1217 0.2 

800818 

124837.7 

1897.5 

100 

18.97 

IS 

800818 

153015. 1 

890818 

161930. 1 

2955.0 

170 

17.38 

10 

009818 

191822.3 

800818 

194952.3 

1890.0 

106 

17.83 

17 

600819 

1 10013.2 

800819 

1 13215. 1 

1919.9 

90 

21.33 

18 

800820 

123830.2 

800820 

133437.7 

2167.5 

44 

49.26 

19 

68G82 1 

1507 7.7 

800821 

154515.2 

2287 . 5 

79 

28.96 

20 

BC0821 

183C30. 1 

80082 1 

191139.9 

2489 . 8 

178 

13.99 

21 

800822 

102622.7 

800822 

105115.2 

1492.5 

38 

39.28 

22 

809322 

134232.8 

800822 

141745.3 

2092.5 

86 

23.78 

23 

809322 

17C915. 1 

800822 

175622.4 

2827 . 3 

172 

16.44 

24 

800323 

92930.2 

800823 

940 0. 1 

629.9 

73 

8.63 

25 

B0C&25 

131730.2 

800825 

134430. 1 

1619.9 

67 

24. 18 

20 

800823 

163237.7 

899825 

.65937.4 

1619.7 

96 

16.37 

27 

800826 

84132.6 

800826 

91352.5 

1919.9 

163 

11.78 

28 

800826 

115322.7 

800826 

122837.7 

1995.0 

191 

10.44 

29 

300826 

153922.7 

800826 

155922.5 

2999.8 

262 

11.45 

30 

300828 

92537.5 

890828 

958 0.2 

1942.7 

97 

20.03 

31 

890828 

1233 0.0 

890828 

131415. 1 

2175. 1 

61 

35.66 

32 

S0C828 

153722.7 

800828 

162230. 1 

1507.4 

41 

36.77 
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Table 36 Observation Summary for Range Data from Station 7090 


Sta t Ion 

« 7090 

Passe 

IS Trackad 

: 30 

Observations total 

4143 

ssassassi 






Baa a aKX a ama a a a ax a 

Pass 

Bearlanlnir Data 

End Ins 
YYMMDD 

Da te 

Pass 

Obs. 

Dens 1 ty 

No. 

YY?IMDD 

Hunnss.s 

Hramss.s 

Dura t Ion 


Lag 






(s) 


(s) 

1 

800614 

73524.0 

800814 

81120.0 

2156.0 

97 

22.23 

2 

800814 

110144.0 

800814 

1130 6.6 

2182.0 

167 

13.07 

3 

800814 

1735 8.0 

800814 

1817 7.0 

2519.0 

182 

13.84 

4 

800814 

210657.0 

800814 

2148 8.0 

2471.0 

207 

11.94 

3 

800815 

93831.0 

800815 

102222.0 

2631.0 

196 

13.42 

6 

800813 

162233.0 

800815 

165245.0 

1810.0 

141 

12.84 

7 

830815 

1941 3.0 

800815 

202838.0 

2853.0 

263 

10.85 

a 

600818 

91445.0 

800818 

94842.0 

2037.0 

1 19 

17. 12 

9 

600818 

124141.0 

800816 

1259 1.0 

1040.0 

67 

15.52 

10 

C00818 

190634.0 

300818 

1953 1.0 

2787.0 

171 

16.30 

1 1 

800819 

111532.0 

800819 

1 14920.0 

2036.0 

136 

14.97 

12 

800019 

175013.0 

800819 

103136.0 

248 1 . 0 

203 

12.22 

13 

800819 

212421.0 

800819 

213335.0 

574.0 

50 

11.40 

14 

830820 

629 6.0 

800820 

65033.0 

1287.0 

29 

44.38 

13 

800820 

2U0224 . 0 

800820 

204123.0 

23.39.0 

136 

17.26 

16 

803821 

83132.0 

000821 

91537.0 

2645 . 0 

104 

25.43 

17 

830822 

72745 . 0 

000322 

74943.0 

1318.0 

55 

23.96 

16 

800822 

104041.0 

800822 

111830.0 

2277 . 0 

162 

14.06 

19 

800322 

171526.0 

80G822 

175657.0 

249 1 . 0 

173 

14.40 

20 

600822 

204736.0 

000822 

2128 5.0 

2409 . 0 

136 

17.71 

21 

600826 

84746 . 0 

800826 

928 7.C 

2421.0 

155 

15.62 

22 

800826 

122016.0 

800826 

123826.0 

1090.0 

41 

26.59 

23 

800826 

1534 1.0 

80C826 

155759.0 

1438.0 

88 

16.34 

24 

800326 

1847 5.0 

300826 

193427.0 

2842.0 

233 

12.20 

23 

800827 

105632.0 

800827 

112919,0 

1947.0 

1 15 

16.93 

26 

SC0827 

143042.0 

BC0827 

143215.0 

93.0 

6 

15.50 

27 

803827 

172746.0 

800827 

101144.0 

2638.0 

189 

13.96 

28 

800827 

210429.0 

800827 

213938.0 

2109.0 

154 

13.69 

29 

630828 

161247.0 

800828 

164735.0 

2088.0 

154 

13.56 

30 

800323 

1937 5.0 

800820 

2022 5.0 

2700 . 0 

214 

12.62 


Both data sets were adjusted using the 6E0SPP program to obtain 
station positions with respect to a fixed orbit. The results for the 
station positions and the associated baseline length from the range and 
SRD adjustments a»e given in Tables 37 and 38 respectively. The two 
baseline estimates differ by about one meter, which considering the fact 
that the SRD observations are good to about 0.5 m and, taking into 
account the sparseness or the data used in this experiment, can hardly 
be used as a basis for drawing firm conclusions about the absolute quali- 
ty of the two estimates. 
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Table 37 Station Coordinates and Baseline Length Estimates and 
Statistics Obtained from the Range Data Adjustment 


Station .'lo . X ( n) Y ( m) Z (m) GX(m) CY( m) GZ(m) 


HL08U2 7090 -2389000.233 3043334.860 -3078527.369 0.043 0.030 0.027 

OHRLAS 7943 -4447344.467 2677140.178 -3694997.100 0.048 0.032 0.033 

Baaellne Leuf^th ( m) 3196328.569 

Standard Deviation ( m) 0.062 


Table 38 Station Coordinates and Baseline Length Estimates and 
Statistics Obtained from the SRD Data Adjustment 


Sta t ion 

No. 

X (n) 

Y (m) 

Z ( n) 

0X(m) 

OY(ni) 

0Z( m) 

nL05O2 

7090 

-2389002.227 

5043333.760 

-3078530.595 

0.744 

0.560 

0.669 

ORRLAS 

7943 

-4447345.783 

2677142.099 

-3693000. 143 

0.595 

0.750 

0.638 

Baae 1 Ine 

Length 

( m) 

3196327.380 





Standard 

Devla t ion ( m) 

0.342 






Besides that, if we consider the location of the stations on the 
earth and the fact that Lageos has a nominal inclination of 109°, we 
reach the conclusion that optimal passes parallel to the dominantly East- 
West direction of this baseline will be hard to come by for this satel- 
lite at any time. We plotted the coobserved events in Fig. 30, and as 
expected almost the entire set of points come from satellite passes 
orthogonal or nearly so to the baseline direction. The deficiency of 
the strongly geometry-dependent SRD mode in such a situation has already 
been pointed out, and it has also been confirmed through the very first 
simulation studies discussed in Section 5.1.1. The results in that case 
(Table 17) indicated that the error of recovery for the ranges would 
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Fig. 30 Lageos SRD Event Distribution for the Data Used in 
Determining the 7090-7943 Baseline 
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only differ (be worse) by 7 mm from that for the SRD solution. If we 
compare this difference to the corresponding differences for some base- 
lines with optimal location with respect to the tracked satellite passes 
(e.g., 7120-7943, 5.6 cm; 7092-7943, 3.2 cm), we realize that it is 
really of an insignificant level. It is not surprising either that the 
results of the second simulation study discussed in Section 5. 1.2.1 also 
point out that this particular baseline is not the best for testing the 
proposed method. Even though in that case the recovery errors for the 
range solution are 3-4 m larger than the corresponding ones for the SRD 
mode, in the case of tue station pair 7090-7943 the difference between 
the two errors is only 18 cm!--hardly significant in the presence of 
10 cm noise. 

The residuals' summary for each of the adjustment are given in 
Tables 39 and 40 for the range solution, and in Table 41 for the SRD's. 
Comparing the mean residual per pass between the two adjustments, we 
find that the SRD solution tends to fit the orbit better for passes 1, 3, 
5, 6, 12, 13 and 14, which, as can be seen from the groundtrack plot in 
Fig. 30, are the ones better conforming with the optimality criteria for 
spatial data distribution in the SRD mode. Finally, a comparison of the 
rms residual between the two solutions shows that the SRD solution tends 
to have residuals with a dispersion which is dictated by the most "noisy" 
of the two stations collecting the observables. In the present case, we 
have already seen in the preliminary adjustment of the complete range 
data set that station 7943 has a noise level 5-10 times higher than sta- 
tion 7090. It is understandable then that the rms residuals in the SRD 
solution are almost identical to those obtained in the range adjustment 
for the data from station 7943. 
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Table 39 Pass-by-Pass Residual Summary for Adjusted Range Data 
from Station 7090 


Pumm 

Ifo. 

Oba. 

Rea Idval 
Mean 
(a) 

rae 

(a) 

Standard 
Devlat ioa 
( a) 

Paaa 

Dorat loa 
(a) 

Hlalaua 
Rea Idual 
(a) 

Hazlana 
Rea Idoa 1 
( a) 

Hean 
Cloaure 
( a) 

1 

83 

- e .«237 

0.088 

0.083 

1049.92 

- 0. 117 

0.097 

- 0.02 

2 

117 

e . eo 38 

0.113 

0.113 

1687.67 

- 0.214 

0.308 

0.00 

3 

37 

- e .2763 

0.291 

0.094 

1064.96 

- 0.418 

- 0.078 

- 0.28 

4 

82 

1879 

0.202 

0.074 

1462.28 

- 0.269 

- 0.039 

- 0 . 19 

3 

55 

e . e ? 9 e 

0 . 128 

0.097 

1192.00 

- 0. 178 

0.238 

0.08 

6 

22 

-o.eess 

0.073 

0.078 

399.00 

- 0. 132 

0 . 146 

- 0.01 

7 

les 

- 0.0361 

1.136 

1.140 

1818.00 

- 1.384 

7.648 

- 0.04 

8 

63 

- 0.0788 

0.099 

0.060 

982.60 

- 0 . 173 

0.059 

- 0.08 

9 

24 

- 0.0120 

0.038 

0.037 

892.60 

- 0. 109 

0.036 

- 0.01 

le 

134 

0.1124 

1.426 

1.426 

2340.04 

- 2.939 

4.611 

0.11 

11 

121 

0. 1202 

0 . 160 

0 . 106 

1884.00 

- 0. 128 

0.384 

0 . 12 

12 

15 

0.0835 

0. 103 

0.068 

473.00 

- 0.082 

0 . 184 

0.08 

13 

87 

0. 1704 

0.212 

0 . 127 

1414.00 

- 0.203 

0.405 

0 . 17 

14 

18 

- 0.4927 

0.497 

0.069 

498.05 

- 0.666 

- 0.440 

- 0.49 


Table 40 Pass-by-Pass Residual Summary for Adjusted Range Data 
from Station 7943 


ana 

No. 

Oba. 

Realdnal 
Hean 
( a) 

raa 

( a) 

Standard 
Devlat Ion 
( a) 

Paaa 

Dora t Ion 
(a) 

Hlniaua 
Rea Idual 
( a) 

HaxlauB 
Rea Idua 1 
( a) 

Hean 
Cloaure 
( a) 

1 

53 

-0.0133 

0.437 

0.441 

1049.92 

- 1 . 040 

0.952 

-0.01 

2 

117 

0.0868 

0.315 

0.304 

1637.67 

-0.716 

0.831 

0.09 

3 

57 

0.0099 

0.393 

0.398 

1064.96 

-0.973 

0.687 

0.01 

4 

82 

0.1142 

0.440 

0.428 

1462.25 

-1. 194 

1.277 

0. 1 1 

3 

53 

0.0634 

0.082 

0.032 

1192.00 

-0 . 043 

0.226 

0.06 

6 

22 

0. 1807 

0.244 

0. 168 

399 . 00 

-0.116 

0.536 

0. 18 

7 

103 

-0.0510 

0.343 

0.341 

1815.00 

-1.098 

0.733 

-0.05 

8 

65 

0.0294 

0.261 

0.261 

932.60 

-0.642 

0.704 

0.03 

9 

24 

-0.0566 

0.578 

0.588 

592.60 

- 1 . 002 

1.449 

-0.06 

10 

154 

-0.0373 

0.260 

0.238 

2340.04 

-0.678 

0.611 

-0.04 

1 1 

121 

-0. 1588 

0. 187 

0.099 

1554.00 

-0.316 

0. 188 

-0. 16 

12 

13 

0.0586 

0.328 

0.334 

473.00 

-1.078 

0.339 

0.06 

13 

87 

0. 1719 

0.203 

0.109 - 

1414.00 

-0.073 

0.360 

0. 17 

14 

18 

-0.0776 

0.453 

0.460 

495.05 

- 1 . 342 

0.574 

-0.08 


oniGlMAL PAGE 13 
OF POOR QUALITY 


180 


ORIGINAL PAGE IS 
OF POOR QUALITY 

Table 41 Pass-by-Pass Residual Summary for Adjusted SRD Data 
from the 7090-7943 Station Pair 


Pass 

No. 

Obs. 

Res idual 
Mean 
( m) 

rms 
( m) 

Standard 
Devlat Ion 
( m) 

Pass 

Dnrat ion 
(s) 

Hlnlanun 
Res Idua 1 
( ■) 

Maximum 
Res Idua 1 
( m) 

Moan 
Closure 
( m) 

1 

53 

-0.0138 

0.414 

0.418 

1049.92 

-1.002 

0.980 

-0.01 

2 

117 

-0. 1076 

0.313 

0.295 

1637.67 

-0.691 

0.591 

-0. 11 

3 

57 

-0. 1331 

0.426 

0.401 

1064.96 

-0.995 

0.694 

-0. 15 

4 

82 

-0.2073 

0.464 

0.418 

1462.23 

- 1 . 363 

1.047 

-0.21 

3 

35 

-0.0381 

0. 176 

0. 173 

1192.00 

-0.381 

0.310 

-0.04 

6 

22 

0. 1796 

0.253 

0. 182 

399 . 00 

-0. 171 

0.480 

0. 18 

7 

103 

-0.0663 

1. 189 

1.193 

1813.00 

-1.680 

7.513 

-0.07 

8 

63 

-0.0985 

0.283 

0.267 

952.60 

-0.714 

0.545 

-0. 10 

9 

24 

0.0382 

0.573 

0.584 

592.60 

- 1 . 467 

0.941 

0.04 

10 

154 

0. 1218 

1.420 

1.419 

2340.04 

-3.544 

4. 152 

0. 12 

1 1 

121 

-0 . 2444 

0.276 

0. 129 

1554.00 

-0.583 

0. 123 

-0.24 

12 

13 

-0.0225 

0.311 

0.321 

473.00 

-1. 130 

0.288 

-0.02 

13 

87 

-0.0774 

0. 134 

0. 110 

1414.00 

-0.364 

0.226 

-0.08 

14 

18 

-0.3999 

0.600 

0.460 

495.05 

-1.010 

0.905 

-0.40 
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6. CONCLUSIONS AND RECOMMENDATIONS 


6.1 Conclusions 

The theoretical investigations and the numerical examples presented 
in this study lead to a number of conclusions concerning the issues 
raised herein. The most important of these is the fact that the analysis 
of Lageos data in the SRD mode minimizes the effect of all of the consid- 
ered types of systematic orbital errors on the estimated baselines and 
coordinates of the pole. We have refrained from attributing these errors 
to any particular source; it is, however, important that we discuss one 
of them. 

It is well known that baseline lengths are independent of the under- 
lying reference system; baseline "estimates," however, especially when 
determined by satellite techniques or even more generally by any non- 
direct measuring system, are directly dependent on the reference system. 

To be more specific, they depend on the stability with which this system 
ca- be maintained. This is a consequence of the fact that the estimate 
is obtained from the end-station coordinates which are determined on the 
basis of their individual observing records. If the "barycenters" of 
these data are considerably apart in time, then the station coordinate 
estimates are affected by the reference system errors accumulated in the 
intervening time interval. The along-track and across-track errors used 
in the simulation study of Section 5. 1.2.1 can be analyzed in a 
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latitudinal and a longitudinal component, and therefore one could think 
of them as the errors in the coordinates of the pole used to rotate the 
satellite positions from the CIS to the CTS frame. 

It is self-evident that when the SRD mode is invoked, the above 
problem is alleviated completely by virtue of the simultaneity of our 
observations. Even though the differencing of the simultaneous ranges 
is not required in this case (simply simultaneous data would be suffi- 
cient), we would recommend that in view of the significant improvements 
in the accuracy of the observable, the SRD mode be followed through in 
its entirety. 

The quality of the results determined on the basis of the recovery 
errors shows that using this method the goal of determining baseline 
lengths with centimeter level accuracies and two-day averages of the 
coordinates of the pole to five centimeters is feasible even in the 
presence of over one meter biases in the orbital model. Such accuracies 
in the orbit are about two to three times our current capabilities in 
predicting the orbit of Lageos over thirty-day periods. 

It is thus conceivable that the predicted Lageos orbit used in 
determining the observational schedules of the stations and the editing 
of the observations can also be used for the analysis of SRD data on a 
nearly real-time basis. The elimination of the satellite orbit from the 
parameter list simplifies the estimation process beyond expectation. 

Users with no access to global sets of data can still use their regional 
data sets in the SPD mode and suffer no loss of accuracy in their results 
even when their reference orbit model is incomplete or they use a fixed 
predicted orbit. The simpler computational procedures of this type of 
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analysis relaxes the expensive hardware requirements imposed by the more 
complicated softwares available. 

The generation of the SRD data set can easily be incorporated into 
the data editing package of the data processing center. As we have seen 
from the results of tests performed on real data, the use of cubic spline 
interpolation is far better than the commonly (and presently) used least 
squares polynomial approximation. Since the latter has several advan- 
tages from the data editing point of view (which is of no concern in this 
study), a compromise between the two methods is probably the best solu- 
tion. Data smoothing cubic splines exist [SpSth, 1974], and in this case 
they would be the most suitable to use. It is our firm belief that the 
above scheme of generating SRD data will result in an insignificant in- 
crease of the overhead cost for data editing which will be well worth it 
considering the improvements in the accuracy of the results and the major 
reduction in the cost of analyzing the data. 

It has been shown here through theoretical arguments as well as 
simulation studies that the SRD mode is very much dependent on the geom- 
etry of the station network and the coobserved Lageos passes. For the 
determination of baseline lengths the best results are obtained from data 
taken on passes which are parallel to the baseline direction. We cannot 
always ensure that such requirements are fulfilled, but we should consider 
doing so whenever we have a choice on the baselines to be determined. For 
all practical purposes, the systematic orbital errors propagate into the 
observables (SRD's) in proportion to the baseline distance separating the 
two coobserving stations. Having a rough estimate of the orbital accuracy 
and the lengths of the baselines between the stations, we can determine 
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the pairs for which the level of the systematic errors remaining in the 
observable after the range differencing will be minimum. The generation 
of bias-surface plots similar to those used in Section 4.4.2 can facili- 
tate the planning stage of a campaign in the selection of the optimal 
station locations as well as the ccobserving pairs of stations. 

In the case of polar motion parameter estimation the situation is 
slightly more complicated. This is mainly because of the fact that not 
only does the station location matter in this case, but the distribution 
of the data in time is of concern here too. We discuss first the station 
location issue. 

As shown in Section 4.5.3 we need two nearly perpendicular baselines 
in order to be able to resolve the two components of the motion of the 
pole. Because of the convention in the definition of these two compon- 
ents X and y, it turns out that the optimal locations are near the two 
prime meridians, i.e., X = 0® or 90° or 180° or 270°. Considering the 
continuous operation of these stations, it is worth pointing out that 
great savings can be achieved if we limit ourselves to an L-shaped 
rather than a +-shaped network, thus decreasing the number of stations 
required to only three. The middle station can be paired with both of 
the outside stations for the required perpendicular baseline pair. 

In an L-shaped network near X = 0° or 180°, the N-S pair is sensi- 
tive to the X component, while the E-W pair determines the y component. 
Exactly the opposite is true in the case of a network near X = ±90°. To 
avoid gravity related orbital errors affecting the estimates, it is ad- 
visable to keep observing networks in both the Northern as well as the 
Southern Hemispheres. 
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The second Issue of interest here is the distribution of the data 
in time. The higher the density of the observation record, the higher 
the resolution in the x and y. In that sense, the resolution is only 
bounded by the accuracy of the observations. The precision of the esti- 
mates can be increased by either increasing the data density over the 
averaging interval or by increasing the interval itself. It should be 
kept in mind that the former does not affect the resolution of the 
parameters, while the latter does. It results in its decrease. These 
are issues to be resolved when a final decision is taken on the practical 
determination of the x and y. In any case though, it should be made 
certain that whatever the chosen averaging interval, there will always 
be enough networks with proper satellite observability schedules to col- 
lect enough well -distributed (globally) observations over the entire 
time interval. 

The analysis of real data has not given us grounds for basing 
firmer conclusions, although the agreement between these results and those 
obtained from the simulation studies gives us a higher degree of confi- 
dence in the validity of the latter. The absence of extensive real data 
tests is due to the lack of suitable data. By this we do not mean to 
cast the blame on others, but rather to point out that as it is shown in 
this study, proper scheduling and a genuine effort from the field parties 
would have certainly resulted in a sizable amount of data. 

There is indeed a striking similarity between the SRD and the pure- 
ly geometric mode. However, the SRD mode requires the coobservation of 
the Lageos pass from only two stations, while in the geometric mode data 
from at least four and preferably more stations with strict simultaneity 
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are required. Apart from the weather factor, there should be no other 
excuse for not obtainino SRD data from nearby stations. The probability 
for two such stations coobserving lageos is much higher than for the 
geometric mode (minimal) four-station network. 

6.2 Recommendations 

The results of this investigation have shown that the proposed 
simultanoous range-differencing approach for the analysis of laser rang- 
ing observations to Lageos is an avenue worth pursuing for improvement 
of our geodetic estimates. On the basis of these results we would recom- 
mend that 

(a) An effort be made during one of the upcoming observational 
campaigns (the 1983 MERIT campaign, for instance) to coobserve as many 
Lageos passes as possible. 

(b) Continue the research effort in optimizing the network config- 
urations that will allow uninterrupted monitoring of the motion of the 
pole within the internationally agreed limits. 

(c) Further research is warranted in the direction of SRD data 
generation. As the field instrumentation is upgraded and the stations 
become capable of obtaining more than one observation per second, the 
amount of incoming data will grow out of proportion. It is therefore 
suggested that further improvement and standardization of the SRD genera- 
tion technique is needed. The possibility of integrating this procedure 
in the raw data preprocessing at the data gathering centers should be 
given serious consideration. 


187 



(d) The scheduling of mobile or highly mobile laser ranging equip- 
ment deployment should be re-examined, and If feasible advantage should 
be taken of this method of analysis by collecting suitable observations. 

(e) As the time for the establishment of a new Conventional Terres- 
trial System nears, the role of the proposed method must be reaffirmed 
through further simulation studies and If possible the analysis of real 
data In contributing optimal estimates for the fundamental polyhedron's 
side lengths. 

(f) A study should be Initiated to Investigate the merits of the 
proposed method In determining the variations in the rotation rate of the 
earth. 

(g) Finally, the application of this method to other ranging or 
pseudo-ranging satellite systems such as the GPS is a research topic 
worth pursuing. 
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APPENDIX A 

DERIVATION OF THE VARIATIONAL EQUATION OF STATE FOR THE 
CASE OF TIDAL ACCELERATIONS 


Equation (29) (Section 2.2.5) gives the perturbing acceleration 
on a satellite due to the tidal effects of a perturbing body b. On the 
basis of this equation we will derive here equation (30) which gives the 
contribution of the above acceleration in the variation of the satellite 
position vector. For clarity we repeat equation (29): 

a ® 

Rjn = I k 2 ([1 - 5(u. *11)2]^ + 2(u. •u)u. ) (A.l) 

'“b ^ iRjjl' |R|- ° ° ° 


where 




(A. 2) 


and 


u,. = 


» 'R. 


To obtain the expression for 
terms in (A.l) individually L 


3R 


TD. 


3U 


3R 


[1 - 5(ujj.u)2]u 


1 - 5 





(A. 3) 

, we differentiate each of the 
-•with respect to R: 

(rJ r)M r 

iRfel^lRI^ 



R 2 


'll' 


31? 

31? 


1 - 5 


(rJ R)' ' 


_ 9|R|'^ 

R ^ + 

3R 
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R|® 


[1 - 5 (Uk*u) 2 ] 


iRj^iRr 
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I + [1 - ^ 


+ < -5 — K 

iRr\ iRbHRr 


2 (Rb^R) <^bL^/,jL 


— -2 

^ V iRr 


[1 - 5(u.-u)2] u J ( )2 _ 

^ I - 5 + 25 —z: u u' 

IRI^ IRI^ IRI^ 

- 10 — ^ u u. ' + 10 _— u u 

\W ^ 1 R|^ 


{[1 - 5{Uh-u)"] I - 5 u + 35 (Ujj-u)" u 

IRr 

- 10 (lTjj'u) u uJ) (A. 4) 


RbTR 


3R \\m/ iRbi 


Vb' 

I^IMRl 


(Rb^R) _ 

+ 2 -IT — "b 

IRbl^ 


-5 1R| — 


= {2(Ub Ub^) - lO(Ub'u) Ub u^} (R-5) 


Collecting terms in the combination of (A. 4) and (A. 5), we obtain: 
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(A. 4) + (A. 5) 


1 


= {[1 - 5(u. -u)2] I + [35(u. -u)2 - 5] u u’’’ 

|R|5 ^ ^ 

- 10(Uj^*u) [uuJ + Uj^u^] + 2(u^J)} 

which when multiplied by the constant terms in (A.l): 

^ "b . 5 


- k 


2 1^1 ’ * 

results in expression (30). 
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APPENDIX B 


SYSTEMATIC CORRECTIONS APPLIED TO THE OBSERVATIONS 

B.l Correction for Geometrodynamical Effects. 

The observed and interpolated ranges are required to compute the 
correction to the SRD's due to signal retardation by gravitation or, 
better, by the curvature of the space [Shapiro et al., 1971; Shapiro, 
1980]. Since the satellite position is also required to compute this 
correction, it seems logical to defer its computation until this position 
is automatically available during the DOC step. 

The determination of the geometrodynamical correction is based on 
the formula that relates TDB end TAI (2.) as obtained by Moyer [1981b] 
and the retardation correction ^or light signals as given in [Shapiro, 
1980]. 

Moyer's expression for aT^ = [TDB - TAI(jl)] can be used to convert 
the measured time interval at the station from a proper time interval to 
a coordinate one, provided we know the epochs that are associated with 
the transmission, reflection and reception of the laser pulse. Since 
though in the case of laser ranging to artificial satellites the whole 
interval rarely exceeds 80 ms, the change of the correction over this 
short time can be obtained from differentiation of AT^ and retention of 
only those terms which are significant. From Moyer's expression for AT^ 

we find that only the second and fourth terms are significant. For a 
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station with cylindrical coordinates (u, X, z), u and z in kilometers, 
we can then write: 

. 1.658 . 10-’ cos E f . 

+ 3.17679 X 10*” u cos(UTl + X) (B.l) 

with 

^ = (1 + e cos M)n (B.2) 

and 


d UTl ^ 
dt “ 


0) 


(B.3) 


where E, M and n are the eccentric anomaly, mean anomaly and the mean 
motion of the earth-moon barycenter's heliocentric orbit respectively, 
and u) is the spin rate of the earth. Substitution in (B.l) results in 


n cos E (1 + e cos M) 

+ 3.17679 X 10’“ 0 ) u cos (UTl + X) (B.4) 

From [Moyer, 1981b] 

n = ^ s 1.99096871 x lO"^ rad/s (B.5) 

and using the adopted value of w for GRS80 [Moritz, 1980] 

u) = 7192 115 X 10"“ rad/s, 
equation (B.4) yields 

“ 3.3010261 X 10"^° cos E(1 + e cos M) 

+ 2.3165518 X 10"“ u cos (UTl + X) (B.6) 

where we have changed the units of u to meters, the final result given 
in seconds. This equation can now be used for the time interval conver- 
sion using the finite observed time interval 6t: 


200 



ORIGINAL PAGE FS 
OF POOR QUALITY 


j(«t) = « (I 
or upon multiplication of the above with c, the speed of light, the 
equivalent correction for the range is 



d(TDB - TAI(Jl)) 
dt 


c6t 


(B.8) 


Since c6t = po. the range prior to the correction, we finally obtain 
6p^ 3 [3.3010261 X 10"“ cos E (1 + e cos M) 

+ 2.3165518 X 10"“ u cos (UTl + X)] po (B.9) 


The remaining retardation correction is computed from 



where and Rj^ are the solar system barycentric coordinates of the 
earth and the satellite respectively. The final range now can be com- 
puted as 

p^ = Po + 6p^ - 6p^ (B.ll) 

the last term being subtracted since it compensates for a retardation. 
With each range in a SRD pair corrected, we can now determine the SRD by 
their difference 

6p = p - p (B.12) 

C L 2 C X 

or using (B.ll) 

6p- = Po 2 + 6p - 6p. - poi - 6p + 6p^ (B.13) 

c rz C 2 ri 

It can ’>e observed though from (B.IO) that 6p^ is independent of 

the station position, and it is the same for both ranges in the SRD pair 

it therefore cancels in the computation of 6p . Some further savings 

c 

can be achieved from a close examination of the 6p^ terms also. From 
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(B.9) 1t can be seen that the second term depends on the station position 
through u, X and po, while the first only on po. With soiiie reasonable 
values for these quantities, the maximal value that 6p^ can reach In the 
case Lageos Is only about ±0.007 m. Upon differencing, therefore. In 
(B.13), the maximum correction for the SRD 6p Is at most ±0.014m, below 
the measuring accuracy of most available Instruments In the field. 

Since the correction hcs a periodic nature [Moyer, 1981b], we can safely 
eliminate It from the computation of 6p^. Equation (B.13) therefore 
takes the simple form 

opj. ~ Po2 “ Poi (B.14) 

which Is used for the determination of the simultaneous range differences 
In this investigation. 

B.2 Systematic Corrections Due to Tidal Motions of the Observing Stations 

The effects of the lunisolar tides on points located on the surface 
of the earth are theoretically rather well studied [Melchior, 1978]. If 
these temporal variations in the location of the observing stations with 
respect to the center of mass of the earth are not accounted for in the 
observations, the committed error can reach an amplitude of about 0.5 m. 

The traditional way of correcting for these effects Is to compute 
the temporal changes in the coordinates of the observing station rather 
than the effect on the measurement directly. In the present study only 
the effect of the solid earth tides was considered mainly due to the 
fact that the remaining effects of the ocean and atmospheric tides are 
much smaller and not yet as well understood or modeled [Lambeck, 1980]. 

The nonrigidity of the earth is accounted for by the Love number for 
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radial expansion h 2 and the Shida number Iz for horizontal shear. The 
values used In our experiments are those used by the NASA/6SFC scientists 
[Chin et al., 1972]: 

hj = 0.600 

I 2 = 0.075 

The mathematical formulation for the station coordinate corrections 
based on the tidal potential from equation (27) Is derived In [Diamante 
and Williamson, 1972]. The accuracy of this formulation can hardly 
match the observation accuracy level today, and it Is well known that 
for best results, the local tidal motions should be obtained from direct 
In situ observation rather than from the model. This has not been the 
case so far for almost all operational SLR stations. 
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APPENDIX C 

FURTHER DEVELOPMENT OF THE DIFFERENTIAL ERROR EQUATIONS 
C.l Station Coordinates 

From the discussion in Section 4.2 of the range differential 
we can use equations (99), (100), and (101) to write 

with 

’'^1/ ■ ^1^ • • V 

Substituting (C.2) in (C.l) we obtain 

= - Ls J [cef - u J [C0][C0]^ ] du. (C. 3 ) 

Pij 1 j j 

From the well-known property of orthogonal matrices R^ = R‘^ 
[Mueller, 1969], we have 

[CO][C0]^ = ,I3 (C.4) 

and (C.3) becomes therefore 

“‘=ij ' ir, - “j"’ 

* w 

Considering now that by (93), 

5/ ■ r/ [Npf (C. 6 ) 

we finally obtain 
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dPij = - [{C0 NP) R. - dUj (C.7) 

^ J 

Note that the first vector in the brackets is the satellite posi- 
tion vector in the CIS frame. The differential relationship for the SRD 
observable d6p^. from stations j,k is obtained by substitution of (C.7) in 
(108). 

C.2 Satellite State-Vector 

The set of equations (99), (100), (102), (103) and (104) yields 
the following range differential relationship: 

‘‘pii = 

Pij LrJ 

Using some of the derivations given in (C.l) and (C.8), we can 

write: 


-[t] 

^ [[NP] R, - [C0]^ Ujf [NPV] d[^-|^ j 


dp,i = [s J - U J [C0]] [NP] [Y 
^ij 


(C.9) 


Using again the property of orthogonal matrices [Ce] [CG] = I, we 
can insert this product between the first and second bracket in (C.9) 
which, upon multiplication, results in 


dp,, = ■— [(C0 NP) R, - U,f [(CC NP) Y] d '1“ 

IJ j I J Ko 


(C.IO) 


The second bracket in (C.IO) is the transitional matrix in the CTS 
frame. The corresponding differential d6p^ for the SRD observations is 
obtained again by differencing (C.IO) written for each of the observing 
stations. 
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C.3 Polar Motion Parameters 

Equations (99), (100), (105), (106) and (107) can be used to write 
the differential relationship between the observed range and the polar 
motion parameters x and y 


* 


’-Uj 3 cose 

-Uj 3 sine 


[C0] 

-Uj 3 sine 

Uj 3 cose 




. ‘'ji 

-“J2 ■ 

(C.ll) 


where we have used some of the substitutions derived in (C.l) and (C.2). 
Multiplying the second with the third matrix we find 



•-Uj 3 cose 

-Uj 3 sine- 

[C0] 

-Uj 3 sine 

Uj 3 cose 


- “jl 

-"j2 - 


-Uj 3 (cos ^0 + sin^Q) + x 
Uj 3 (cos 0 sin6 - cosQ sin6) - y Ujj 


L Uj 3 (x cos^Q + y sin0 cos9 + x sin^0 - y sin0 cos0) + U 


jl 


-Uj 3 (cos 0 sin0 - sin0 cos0) 


Uj3(sin"e + cos^e) 


Uj 3 (x cose 

sine + y 

CD 

CM 

C 

j3 

- 

* V 


'J3^ 


jl " * «j3 


>^3. 


The last approximation in (C.12) 
lation studies, even in actual solutio 


X sin0 COS0 + y cos^0) - 0^2 



■-'^03 

0 ■ 

= 

0 

“j3 


.“jl 

— — 1 

CVJ 

ZD 

1 


is justified in the case of simu- 
s indeed, since the x,y parameters 
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are of the order of 10"® radians and the stations are confined on the 
earth's surface. Therefore, 

|UjJ - 0 (10®) (C.13) 

To obtain the differential d6p^. for the SRD observation, we must 
evaluate (C.ll) for the two coobserving stations and subtract the result- 
ing expressions. 
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RESIDUAL SUMMARIES FOR TEN SLR STATIONS FOR 
THE AUGUST^ 1980, LAGEOS DATA 
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Table 42 Residual Sunmary for Station 7063 


CONSOLIDATED STATISTICS FO« STATION I TOAl 


FASS 

08SERV 

aesio MEAN 

RNS 

DEVIATION 

length 

NIN AESO 

NAX AESO 

MEAN CLCS 

1 

4 

B WB A BBB B B BBBBBB i 

-0.«2tt 

S.MS 

BB BBS BBB B BBI 

4.123 

1292.00 

-8.493 

3.007 

-0.93 

2 

477 

0.1070 

0.237 

0.212 

1441.00 

-2.866 

0.398 

0.11 

i 

202 

0.1237 

0.443 

0.432 

1494.00 

-4.840 

3.480 

0.12 

A 

4 

>2.0974 

3.439 

3.012 

1689.00 

-3.903 

1.408 

-2.10 

S 

• S9 

0.1242 

0.223 

0.187 

2338.00 

-2.436 

0.473 

0.12 

A 

1 

0.04S4 

0.046 

0.0 

0.0 

0.046 

0.046 

0.05 

T 

ISSO 

0.0119 

0.322 

0.321 

2810.00 

-4.383 

8.987 

0.01 

• 

4 

-4.4022 

3.623 

4.043 

1303.00 

-9.652 

-1.043 

0 

1 

9 

14 

-0.4982 

2.473 

2.314 

2330.00 

-4.545 

3.946 

-3.50 

ID 

1167 

-0.1706 

0.444 

0.432 

2484.00 

-6.694 

7.124 

-C.17 
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CONSOLIDATED STATISTICS 

FOA station 

t 7090 





FASS OeSEAV AESID NEAS 

BBBBBBBBBBBBBBBBBBr •••••BBBBBBBBB 

RNS 

BBBBBB BBBBBB 

deviation 

length 

MIN MESD MAH RESO MEAN CLOS 

BBBBBLBBBBBBBBBBBBBBBBBBBBBBBBBB 

1 

97 

0.0882 

0.130 

0.095 

2136.00 

-0.221 

0.144 

0.09 

2 

167 

-0.0325 

0.104 

0.099 

2182.00 

-0.117 

0.177 

-0.01 

3 

182 

-0.0892 

0.131 

0.096 

2319.00 

-0.513 

0.111 

-0.09 

4 

207 

0.0322 

0.109 

0.105 

2471.00 

-0.282 

0.264 

0.01 

5 

196 

0.0135 

0.140 

0.117 

2631.00 

-0.431 

0.44 8 

0.03 

4 

141 

-0.0832 

0.142 

0.1 *6 

1810.00 

-0.386 

0.192 

-0.08 

7 

263 

-0.0764 

0.119 

0.091 

2851.00 

-0.427 

0.177 

-0.08 

8 

119 

0.0891 

0.118 

0.106 

2017.01 

-0.151 

0.384 

0.09 

9 

67 

-0.0531 

0.093 

0.078 

1040.00 

-0.215 

0.145 

-0.05 

10 

i;i 

-0.0780 

0.122 

0.094 

2787.00 

-0.471 

0.181 

-0.08 

II 

116 

-0.0981 

0.118 

0.097 

2016.00 

-0.451 

0.101 

-0.10 

12 

203 

-0.0640 

0.112 

0.092 

2481.00 

-0.526 

0.145 

-0.36 

13 

50 

0.0121 

0.080 

0.079 

5 74.00 

-0.157 

0.216 

0.01 

14 

29 

0.0674 

0.141 

0.128 

1287.00 

-0.219 

0.286 

0.07 

15 

116 

-0.0940 

0.124 

0.081 

2319.00 

-0.147 

0.092 

-0.09 

16 

104 

0.1106 

0.415 

0.422 

2645.00 

-1.875 

0.378 

0.11 

IT 

55 

0.1305 

0.213 

0.194 

1118.01 

-0.282 

0.693 

0.13 

18 

162 

-0.1606 

0.204 

0.126 

2277.00 

-0.472 

0.071 

-0.16 

19 

171 

U.0T15 

0.116 

0.114 

2491.00 

-0.272 

0.126 

0.07 

20 

lie 

-0.1598 

0.212 

C.140 

2439.00 

-0.567 

0.111 

-0.16 

21 

155 

0.0941 

0.185 

0.159 

2421.00 

-0.321 

0.178 

0.09 

22 


0.0111 

0.069 

0.062 

1090.00 

-0.111 

0.141 

0.03 

21 

88 

0.1992 

0.218 

0.088 

1438.00 

-0.102 

0.148 

0.20 

24 

211 

-0.0411 

0.112 

0.104 

2842.00 

-0.597 

0.195 

-0.04 

25 

115 

0.1264 

0.156 

0.092 

1947.00 

-0.168 

0.145 

0.13 

26 

6 

0.1555 

0.172 

0.071 

93.00 

0.071 

0.280 

0.16 

27 

189 

0.0197 

0.128 

0.126 

2618.00 

-0.289 

0.294 

0.02 

28 

154 

0.08S7 

0.121 

0.085 

2109.00 

-0.164 

0.128 

0.09 

29 

154 

-0.0249 

0.223 

0.221 

2088.00 

-0.567 

0.198 

-0.02 

30 

214 

0.0598 

0.128 

0.114 

2700.00 

-0.179 

0.292 

0.06 
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Table 44 Residual Sunmary for Station 7091 


CONSOIIOATEO STATISTICS FOR STATIOM i T091 
PASS OeSERV RESIO MEAN RMS DEVIATION 


length min reso max reso mean clos 


1 

m 

0.08ZA 

0.2T1 

0.2S9 

HAS. 00 

-O.TSO 

0.4R] 

0.04 

2 

3S2 

-0.06A4 

0.1S2 

O.ITI 

ITli.OO 

-l.OAO 

0.1A4 

-0.04 

i 

ZAO 

O.OSOS 

0.149 

0.142 

IRT9.00 

-0.9A5 

O.ASO 

0.05 

A 

A)9 

-0.0202 

O.llA 

O.IIS 

194S.0I 

-1.IA2 

A.TA2 

-0.02 
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consol 

lOATED STATISTICS FOR STATION 1 T092 





PASS OBSERV 

HESIO NEAN R*tS DEVlAriON 

length 

MIN RESO 

MAX RCSD 

sama 

MEAN CLOS 

■ ■•••aaaaai 

1 

122 

0.2185 

0.244 0.154 

IT41.01 

-0.258 

0.52T 

0.22 

2 

246 

-0.IT12 

0.248 0.204 

1 1 53. 00 

-0.8A5 

0.444 

-O.IT 

S 

12T1 

-O.OOOA 

0.239 0.239 

2180.00 

-1.21A 

0.984 

-0.00 

4 

143 

-0.032A 

O.lOA 0.101 

22T1.0C 

-1.291 

O.TIA 

-0.01 

5 

9 

0.0924 

0.131 0.13T 

ITA.OO 

-O.MT 

0.T26 

0.09 
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CONSOLIDATED STATISTICS 

FOR STATION 

t T094 





PASS 

OBSERV 

RESIO MEAN 

RMS 

la 

DEVIATION 

LENGTH 

MIN RESO MAX RESO MEAN CLCS 

1 

949 

O.OOT8 

0.189 

0.189 

2189.00 

-0.581 

0.5A4 

0.01 

2 

A4I 

0.0359 

0.150 

0.1A4 

2008.99 

-O.Tll 

0.111 

O.OA 

3 

248 

-0.1355 

0.25T 

0.219 

1109.01 

-0.931 

0.313 

-O.IA 

A 

91 

-O.IOTS 

0.391 

0.2AA 

452.00 

-0.953 

0.122 

-0.11 

5 

A5 

0.05AT 

0.144 

0.158 

92A.00 

-0.A51 

0.354 

0.06 

4 

416 

0.0351 

0.213 

0.210 

1148.01 

-1.019 

0.531 

O.OA 
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CONSOLIDATED STATISTICS 

FOR STATION 

t TllA 





PASS 

08SCRV 

RESIO MEAN 

RMS 

■aaaaaaaaaaaa 

DEVIATION 

aaaaaaaasaa 

LENGTH 

MIN RESO MAX RESO MEAN CLOS 

1 

182 

-0.0A05 

0.IT4 

0.1T2 

1181.99 

-0.A58 

1.001 

-O.OA 

2 

IT 

-0.13A4 

I.AT5 

1.51A 

1114.00 

-A.9T9 

2.412 

-0.11 

1 

855 

0.0155 

0.241 

0.242 

2535.00 

-3.945 

1.944 

0.02 

A 

9 

1.1201 

2.A90 

2.150 

114T.00 

-1.T29 

5.192 

1.12 

5 

141 

0.0919 

0.155 

0.I2A 

1009.99 

-0.110 

0.A45 

0.09 

4 

190 

-0.0014 

0.129 

0.129 

2102.99 

-0.18T 

0.948 

-0.00 

T 

4 

-0.T158 

1.455 

1.922 

88T.00 

-A.818 

4.0T2 

-0.7A 

8 

228 

-0.0292 

0.150 

0.1A9 

10A5.00 

-0.18A 

A. 294 

-0.03 

9 

T 

O.OTTl 

0.111 

0.084 

4T4.00 

0.001 

0.214 

0.08 

10 

T 

0.A2A1 

O.TOl 

0.404 

1080.00 

0.099 

I.TT9 

0.A2 

11 

A 

-A. 5052 

5.T95 

A.208 

AOO.OO 

-8.094 

3.451 

-A. 51 


211 


ORIGINAL PAGE 13 
OF POOR QUALITY 
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CONSOLIDAUO STATISTICS FOR STATION I T119 


FASS 

OBStRV 

RESIO NEAN 

RNS 

DEVIATION 

lENCTH 

NIN RESD 

PU RESD 

NEAN CLOS 

I 

266 

0.OB9O 

O.IBI 

0.099 

1226.01 

-0.178 

0.179 

0.08 

2 

29 

0.2660 

0.280 

e.OBB 

711.00 

0.077 

0.610 

0.27 

1 

166 

-O.OTIO 

1.09T 

1.096 

2001.00 

-6.960 

0.999 

-0.07 

M 

27 

-6.7929 

6.79} 

0.090 

1021.00 

-6.969 

-6.61T 

-6.79 

i 

900 

0.0916 

0.191 

0.119 

2271.00 

-0.622 

1.167 

0.09 

6 

sa 

0.1169 

0.328 

0.08B 

1911.00 

0.120 

0.688 

0.12 

7 

172 

0.1989 

0.911 

0.166 

1166.00 

-0.929 

6.826 

0.16 

S 

61 

-0.0666 

0.169 

0.118 

889.00 

-0.166 

0.297 

-0.05 

9 

119 

-0.2916 

0.112 

O.IBI 

1911.00 

-0.719 

0.906 

-0.25 

10 

988 

0.1 TOB 

0.211 

0.I2B 

2727.00 

-0.299 

1.811 

0.17 

U 

17 

0.1109 

0.169 

0.102 

1608.00 

-0.079 

0.161 

0.11 

12 

44 

0.2616 

0.286 

0.191 

692.00 

-0.060 

0.969 

0.26 

Table 

49 

Residual Summary for Station 

7120 





CONSOLIOATEO STATISTICS 

FOR STATION 

t TI20 





PASS 

OBSERV 

RESIO NEAN 

RNS 

DEyiV!9N_ 

LENGTH 

NIN RESO 

MAX RESO 

NEAN CLOS 

1 

229 

-0.1211 

0.162 

0.071 

1966.00 

-0.180 

0.120 

-0.12 

2 

66 

-O.OOST 

0.197 

0.199 

829.00 

-0.772 

0.180 

-0.00 

S 

160 

0.0996 

0.160 

0.098 

1297.00 

-0.229 

0.111 

0.10 

4 

62 

-0.0669 

0.108 

0.086 

618.00 

-0.218 

0.098 

-0.07 

i 

187 

0.026B 

0.111 

0.111 

2616.00 

-0.117 

0.857 

0.01 

A 

166 

-0.0766 

0.116 

0.089 

2799.00 

-0.168 

0.267 

-0.09 

7 

601 

0.0911 

0.118 

0.102 

2971.00 

-0.299 

0.298 

0.09 

a 

90 

-0.29B1 

0.296 

0.161 

869.00 

-0.911 

0.016 

-0.26 

4 

121 

-0.1 BT9 

0.216 

O.IOT 

1699.00 

-0.180 

0.106 

-0.19 

10 

128 

0.1102 

0.161 

0.120 

2617.00 

-0.221 

0.609 

0.11 

Table 

50 

Residual Summary for Station 

7907 





CON SOI 

lOATEO STATISTICS 

FOR STATION 

1 T90T 





PASS 

OBSERV 

4ES1D UfAN RMS 

DEVIATION 

LENGTH 

NIN RESD 

■•■••■•••a 

NAX RESD 

NEAN CLOS 

aaa«aaaaaai 

1 

12 

0.1211 

0.166 

0.160 

292.60 

-0.608 

0.919 

0.12 

2 

61 

-0.0119 

0.606 

0.611 

1112.99 

-1.920 

0.761 

-0.01 

S 

91 

0.0917 

0.909 

0.909 

862.91 

-1.291 

0.801 

0.09 

4 

19 

-0.0686 

0.119 

0.161 

1027.76 

-0.767 

0.961 

-0.09 

S 

5 

-0.0086 

0.996 

0.666 

160.19 

-1.166 

0.666 

-0.01 

6 

92 

0.1667 

0.921 

0.698 

892.90 

-0.991 

2.012 

0.17 

7 

26 

0.0601 

0.601 

0.609 

660.27 

-0.929 

0.726 

0.06 

a 

19 

- 0 . 202 a 

0.191 

0.166 

607.66 

-0.868 

0.160 

-0.20 

9 

19 

0.0691 

0.288 

0.287 

1162.96 

-0.64B 

0.660 

0.09 

10 

91 

-0.0616 

0.^08 

0.608 

916.99 

-1.008 

0.881 

-0.06 

11 

16 

-0.0226 

o.S9i; 

0.601 

967.91 

-1.199 

0.697 

-0.0? 

12 

5 

0.2901 

0.297 

0.069 

160.10 

0.169 

0.101 

0.29 

IS 

26 

-0.9918 

1.009 

0.810 

1110.02 

-2.061 

0.962 

-0.99 

1« 

IT 

0.0006 

0.811 

0.818 

689.98 

-1.916 

1.969 

0.00 

IS 

B 

0.2717 

0.928 

0.686 

620.00 

-0.616 

1.126 

0.27 

16 

9 

-0.1126 

0.968 

0.990 

669.09 

-1.112 

0.692 

-0.11 

IT 

22 

-0.1899 

0.766 

0.T17 

867.96 

-1.9»0 

1.861 

-0.19 

18 

28 

0.2187 

0.609 

0.111 

1012.90 

-0.646 

0.796 

0.26 

19 

2 

-0.198 7 

0.629 

0.207 

90.00 

-0.969 

-0.292 

-0.60 

20 

29 

0.0272 

0.61 1 

0.617 

779.97 

-0.899 

0.698 

0.01 
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Table 51 Residual Summary for Station 7943 



CONSOL lOATEO 

STAIISTICS 

FOR STATION 

t 7963 





PASS OBSEAV RESIO MEAN 

•■■■■■a 

RMS 

OCVIATION 

LENGTH 

NIN RESO 

aaajtaasaa m 

NAR RESO 

•••■maaaaa 

MEAN CLOS 

1 

96 

0.0012 

0.633 

0.637 

1196.92 

-1.039 

0.999 

0.00 

2 

122 

0.0626 

0.338 

0.336 

2329.29 

-1.632 

0.819 

0.06 

1 

82 

-0.0668 

0.639 

0.637 

1909.10 

-1.292 

1.076 

-0.07 

4 

98 

-0.0301 

0.608 

0.609 

2272.36 

-0.970 

0.719 

-0.03 

s 

83 

0.0633 

0.626 

0.626 

1686.79 

-1.229 

1.222 

0.06 

A 

89 

0.0916 

0.337 

0.339 

1967.79 

-0.917 

0.872 

0.09 

T 

90 

0.1379 

0.292 

0.299 

1626.99 

-0.689 

0.991 

0.16 

• 

212 

-0.0796 

0.318 

0.309 

2766.99 

-0.866 

0.679 

-0.08 

« 

199 

0.0983 

0.366 

0.330 

2910.00 

-0.720 

0.861 

0.10 

10 

67 

-0.0263 

0.911 

0.916 

1936.86 

-1.092 

1.198 

O 

o 

11 

96 

0.0337 

0.629 

0.632 

1916.69 

-1.193 

1.393 

0.03 

12 

79 

-0.1026 

0.992 

0.966 

1609.00 

-1.306 

0.976 

e 

0 

O 

1) 

139 

0.0616 

0.239 

0.236 

1822.68 

-0.681 

0.996 

0.06 

14 

100 

-0.0036 

0.379 

0.377 

1897.69 

-0.896 

1.061 

-0.00 

IS 

170 

-0.0620 

0.616 

0.613 

2999.00 

• 

o 

1.299 

O 

O 

16 

106 

-0.0837 

0.367 

0.338 

1890.00 

-1.137 

0.719 

-0.08 

IT 

90 

-0.0132 

0.280 

0.281 

1919.92 

-0.878 

0.702 

-0.01 

18 

44 

-0.0889 

0.670 

0.667 

2167.67 

-1.129 

0.663 

-0.09 

19 

79 

0.1667 

0.626 

0.392 

228 7.99 

-0.963 

0.899 

0.17 

20 

178 

0.0212 

0.282 

0.282 

2689.80 

-0.607 

1.293 

0.02 

21 

38 

-0.1376 

0.629 

0.617 

1692.99 

-1.638 

1.660 

O 

22 

88 

-0.0669 

0.689 

0.683 

2092.96 

-1.167 

1.269 

-0.07 

23 

172 

-0.0399 

0.263 

0.261 

2827.32 

-0.»71 

0.6 1 

'0.06 

24 

73 

-0.3093 

0.399 

0.299 

629.96 

-0.896 

0.119 

-J.31 

29 

67 

0.0391 

0.668 

0.690 

1619.87 

-1.126 

0.798 

o 

• 

o 

26 

96 

-0.0969 

0.996 

0.992 

1619.69 

-1.609 

1.369 

-0.09 

27 

163 

-0.1619 

0.279 

0.261 

1919.90 

-0.771 

0.911 

-0.16 

28 

191 

0.1109 

0.330 

0.312 

1996.97 

-0.72 7 

0.898 

0.11 

29 

26? 

0.1799 

0.369 

0.318 

2999.78 

-0.819 

0.837 

0.18 

30 

97 

0.1617 

0.282 

0.269 

1962.69 

-0.673 

0.703 

0.16 

31 

61 

0.1398 

0.666 

0.668 

2179.10 

-0.993 

0.879 

0.16 

32 

61 

-0.0869 

0.689 

0.683 

1907.66 

-1.391 

0.607 

-0.09 
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